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Chapter  1 
INTRODUCTION 

1 . 1  Problem  Background 

The  two  standard  computational  techniques  for  treating  boundary- 
layer  flows  are  the  direct  and  inverse  methods.  The  direct  method, 
which  involves  specifying  a  pressure  or  inviscid  edge  velocity 
distribution,  has  been  widely  used  for  many  years.  The  pioneering  work 
of  Smith  and  Clutter  [1]  is  an  excellent  example  of  a  direct  boundary- 
layer  solution  method.  This  approach  performs  well  for  attached  shear 
layers  but  cannot  treat  separating  boundary  layers,  due  to  the 
singularity  that  exists  at  the  separation  point  when  the  edge  velocity 
is  specified  [2,3].  This  singularity  prevents  the  use  of  the  direct 
method  for  regions  of  backflow,  including  separation  bubbles  and 
trailing  edge  separation. 

These  regions  of  backflow  are  usually  of  great  interest  to  the 
fluid  dynamicist,  since  they  can  have  a  significant  effect  on  the 
overall  flow.  To  obtain  solutions  for  separated  boundary  layers, 
several  types  of  inverse  methods  have  been  developed.  By  specifying 
distributions  of  displacement  thickness,  wall  shear,  or  similar 
quantities,  and  obtaining  the  streamwise  pressure  distribution,  these 
methods  eliminate  the  singularity  in  the  boundary- layer  equations  at 
separation.  The  equations  of  motion  may  then  be  integrated  through  the 
region  of  separation  under  the  condition  that  the  extent  of  separation 
is  small  (backflow  velocities  less  than  ten  percent  of  the  edge 
velocity).  The  boundary- layer  equations  remain  valid  under  the 
assumption  that  the  shear  layer  is  thin  (d5/dx  <<  1).  For  shear  layers 
which  violate  this  condition,  upstream  influences  become  significant  and 
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the  boundary- layer  equations  alone  can  no  longer  accurately  describe  the 
flow. 

Catherall  and  Mangier  [4]  originated  the  inverse  procedure, 
specifying  a  displacement  thickness  distribution  near  the  separation 
point  to  obtain  regular  solutions  into  a  region  of  backflow.  Once 
inside  the  separated  zone,  however,  the  procedure  developed 
instabilities  which  required  a  progressive  increase  in  convergence 
limits  as  the  solution  continued  downstream.  Since  the  boundary- layer 
equations  are  parabolic,  these  instabilities  arise  from  marching  in  a 
direction  opposite  to  the  streamwise  velocity. 

Reyhner  and  Flugge-Lotz  [5]  originated  the  FLARE  approximation 
which  eliminates  the  instabilities  by  setting  all  streamwise  convection 
terms  to  zero  inside  the  region  of  backflow.  Since  the  x-convective 
terms  are  small  in  this  region,  for  thin  separated  zones,  this 
introduces  only  minor  errors  while  allowing  a  forward-marching  procedure 
to  continue  through  the  separated  region.  Since  it  was  presented,  the 
FLARE  approximation  has  been  a  standard  part  of  the  majority  of  inverse 
methods  which  seek  to  resolve  regions  of  separated  flow. 

An  iterative  procedure  developed  by  Williams  [6]  improved  on  the 
FLARE  approach  by  making  repeated  downstream  and  upstream  passes  through 
the  separated  zone.  On  the  first  downstream  pass,  FLARE  is  used  to 
define  the  extent  of  the  separation  zone.  An  upstream  pass,  using  a 
direct  solution,  is  then  made  inside  the  region  of  backflow  to  determine 
the  convective  terms.  The  downstream  pass  is  then  repeated  without 
FLARE,  using  the  results  obtained  from  the  upstream  solution.  This 
process  is  repeated  for  several  cycles  until  convergence  is  reached. 

The  DUIT  (Downstream-Upstream  Iteration)  procedure  produces  more 
accurate  velocity  profiles  for  the  separated  region,  and  has  been  used 
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successfully  by  Cebeci,  Keller,  and  Williams  [7].  A  drawback  inherent 
in  the  iterative  process,  however,  is  that  only  closed  separation 
bubbles  or  separated  regions  which  reach  an  asymptotic  state  can  be 
treated.  The  DUIT  procedure  has  presently  been  tested  only  on  laminar 
flows . 

The  majority  of  inverse  studies  have  only  considered  laminar 
boundary  layers.  These  studies  include  the  work  of  Keller  and  Cebeci 
[8] ,  Klineberg  and  Steger  [9],  Horton  [10,11],  and  Carter,  [12,13]. 
Recently,  several  studies  of  turbulent  separating  boundary  layers  have 
been  presented,  including  Carter  and  Womom  [14],  Cebeci  [15],  Carter 
[16],  and  Whitfield,  Swafford,  and  Jacocks  [17].  These  are  difference 
methods,  with  the  exception  of  the  integral  method  in  Reference  17. 

The  quantities  specified  by  these  inverse  methods  vary,  with  most 

specifying  the  Reynolds  stretched  displacement  thickness  [4]  or  a  wall 

shear.  A  method  presented  by  Carter  [16]  is  unique  in  specifying  the 

mass  flow,  m.  Carter  [12,13]  has  considered  the  specification  of  both 

displacement  thickness  and  wall  shear,  but  reformulated  the  equations 

for  each  case.  Cebeci  [15],  however,  considers  the  specification  of 

skin  friction  coefficient  and  displacement  thickness  with  a  minimum 

e 

of  reformulation.  Cebeci  also  presents  both  direct  and  inverse  methods 
in  this  work,  although  the  equations  must  be  put  in  variational  form  to 
obtain  an  inverse  solution. 

Table  1  presents  a  sample  of  the  inverse  methods  referenced  in  this 
study,  with  a  brief  summary  of  their  important  features. 


Table  1 


Sampling  of  Inverse  Method  Investigations 


Authors 

Specify 

Flow  Tvpe 

Method/Discretization 

1 .  Catherall  and 
Mangier  (1966) 

6* 

Lam 

Crank-Nicolson 

2.  Reyhner  and 

Flugge-Lotz  (1968) 

JU 

6~ 

Com/ Lam 

Implicit  F.D. /Integral 

3.  Keller  and 

Cebeci  (1972) 

1  t 

f 

w 

Lam 

Eigenvalue/Keller  Box 

4.  Cebeci  and 

Keller  (1973) 

i  i 

f 

w 

Lam 

Mechul  Fun. /Keller  Box 

5.  Klineberg  and 
Steger  (1974) 

T 

Lam 

Modified  Box 

6.  Horton  (1974) 

*r 

fw 

Lam 

Shooting 

7.  Carter  (1974) 

V** 

Lam 

Modified  Box 

8.  Carter  (1975) 

Vs* 

Lam 

Global  It. /Modified  Box 

9 .  Carter  and 

Wo  mum  (1975) 

<5* 

Lam/Tur 

Modified  Box 

10.  Cebeci  (1976) 

Cf/5* 

Com/Lam/Tur 

Variational  Eq. 

Keller  Box 

11.  Carter  (1978) 

6*/m 

Com/Lam/Tur 

Crank-Nicolson 

12.  Cebeci,  Keller,  6 

and  Williams  (1979) 

Lam 

DUIT/Keller  Box 

13.  Horton  (1981) 

A 

6" 

Lam 

Implicit  F.D. /Imbedding 

14.  Whitfield, 

Swafford,  and 
Jacocks  (1981) 

6* 

Lam/Tur 

Integral 

ivttf 


**  . 
k\ 
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1 . 2  Applications 

Before  discussing  the  present  approach  to  solving  the  boundary- 
layer  equations,  a  brief  consideration  of  the  applications  of  direct  and 
inverse  boundary- layer  methods  is  appropriate.  If  the  solution  method 
cannot  be  applied  to  important,  existing  problems,  any  innovations  in 
the  solution  approach  are  worthless. 

As  is  stated  in  the  previous  section,  the  direct  method  is 
applicable  only  to  attached  flows.  It  is  typically  used  to  obtain  skin 
friction  distributions,  drag  coefficients,  and  integral  parameters  which 
characterize  the  viscous  flow  over  a  body  or  surface  subjected  to  a 
particular  pressure  distribution.  Because  of  the  familiarity  of  the 
direct  method,  we  will  concentrate  this  discussion  on  the  application  of 
inverse  methods. 

Inverse  methods,  which  apply  to  both  attached  and  separated 
boundary  layers,  provide  a  pressure  distribution  when  supplied  with  a 
displacement  thickness  or  wall  shear  distribution.  Inverse  methods  have 
several  important  applications  when  used  alone,  and  can  also  be  used  in 
conjunction  with  direct  boundary- layer  methods  and  inviscid  solution 
methods. 

Used  alone,  inverse  methods  can  be  effective  design  tools.  This  is 
particularly  true  of  the  wall  shear  methods,  where  an  optimal  design  for 
body  shapes  can  be  developed  by  prescribing  the  friction  forces  on  a 
surface  and  determining  the  corresponding  pressure  distribution.  Such 
methods  could  also  be  used  for  evaluating  flow  situations  where  the 
boundary  layer  maintains  near  zero  wall  shear  values.  An  example  of 
this  situation  is  the  Liebeck  airfoils  [18],  which  are  designed  for  high 
lift  characteristics. 
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An  important  application  of  inverse  methods  is  interacting 
boundary- layer  methods  or  viscous/ inviscid  interaction  schemes.  When 

-1/2 

the  boundary- layer  thickness  remains  of  order  (Re)  ,  the  flow 
behavior  downstream  has  negligible  effect  on  the  upstream  flow.  If  the 
boundary- layer  thickness  becomes  order  unity,  however,  the  upstream 
pressure  effect  is  significant  and  the  boundary- layer  equations  no 
longer  adequately  describe  the  flow.  In  these  situations,  the  boundary- 
layer  solution  is  usually  coupled  to  an  inviscid  outer  solution. 

Although  this  approach  is  useful  for  regions  of  weak  interaction,  with 
minor  pressure  effects,  it  is  particularly  appropriate  for  regions  of 
strong  interaction  with  boundary- layer  separation.  This  is  a 
satisfactory  alternative  to  solving  the  full  Navier-Stokes  equations. 

This  area  of  computational  fluid  dynamics  is  too  complex  to  be 
thoroughly  discussed  here,  with  many  variations  in  the  application  of 
boundary- layer  solutions,  both  direct  and  inverse.  For  our  purpose,  it 
is  sufficient  to  make  the  reader  aware  of  the  use  of  boundary- layer 
solutions  in  this  manner.  For  a  more  complete  discussion  refer  of 
Veldman  [19]  or  Vatsa  and  Verdon  [20]. 

1 . 3  The  Present  Method 

This  study  presents  a  unified  approach  to  the  numerical  solution  of 
the  direct  and  inverse  formulations  of  the  boundary- layer  equations 
using  an  implicit  spline/finite  difference  discretization  scheme.  The 
method  is  formulated  to  allow  the  specification  of  transformed  wall 
shear,  skin  friction  coefficients,  or  displacement  thickness 
distributions  by  simply  switching  the  appropriate  boundary  conditions. 
This  idea  is  similar  to  the  work  of  Cebeci  [15].  The  current  study  goes 
beyond  that  of  Cebeci  by  also  allowing  direct  or  inverse  solutions  to  be 


obtained  with  a  similar  technique.  Thus  the  present  work  is  truly 
unified  in  handling  direct  as  well  as  shear  or  displacement  thickness 
inverse  solutions  with  the  same  formulation,  requiring  only  minor 
changes  in  the  boundary  conditions  depending  upon  the  desired  solution. 

The  present  scheme  applies  fourth  order  splines,  as  derived  by 
Rubin  and  Khosla  [21],  to  approximate  the  normal  derivatives  in  the 
boundary-layer  equations,  with  two-and  three-point  backward  differences 
to  discretize  the  streamwise  derivatives.  The  differencing  scheme  was 
chosen  to  yield  a  fully  implicit  solution  method.  The  splines  also 
provide  increased  accuracy  for  a  given  number  of  mesh  points. 

The  mechul  function  method,  originally  devised  for  inverse  methods 
by  Cebeci  and  Keller  [22],  is  used  in  a  modified  form  described  by 
Kaufman  and  Hoffman  [23].  This  modified  form  is  necessary  due  to  the 
use  of  the  spline  discretization.  The  FLARE  approximation  is  employed 
in  separation  zones  to  prevent  the  development  of  instabilities  inherent 
in  marching  downstream  through  regions  of  backflow. 

To  treat  turbulent  boundary  layers,  a  two-piece  algebraic  eddy 
viscosity  model,  developed  by  Baldwin  and  Lomax  [24],  is  incorporated  in 
the  formulation.  The  modified  Levy-Lees  transformation  is  used  to 
capture  the  boundary- layer  growth  for  laminar  and  turbulent  flows.  A 
modified  version  of  the  adaptive  normal  coordinate  developed  by  Carter, 
Edwards,  and  Werle  [25]  is  used  to  resolve  the  wall  region  of  the 
turbulent  boundary  layer. 

Following  discretization,  the  nonlinear  spline/finite  difference 
equations  are  linearized  and  solved  using  Newton's  method.  The 
resulting  linear  block  tridiagonal  matrix  system  for  the  solution  vector 
corrections  at  each  streamwise  station  is  solved  using  L-U 


decomposition.  Partial  pivoting  is  necessary  in  the  block  tridiagonal 
solution  process  to  prevent  the  buildup  of  roundoff  errors. 


Chapter  2 
ANALYSIS 
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2. 1  Governing  Equations 

The  governing  equations  for  steady,  incompressible,  two- 
dimensional  boundary- layers  in  dimensionless,  Reynolds  number  scaled 


form  are 
continuity 

x-momentum 

y-momentum 

where. 


*"  +  *1  =  0 
ax  3Y 


3u  3u 

PU  Ti  +  PV  3Y  ' 


3p  ,  3  i  3un 
3x  +  3Y  (yt  3Y 


=  0 
3Y 


=  y  +  pe 


and  the  Euler  equation  applied  at  the  body  surface  gives 

dp  du 

*e  _  _ e 

dx  pe  Ue  dx 


(1) 

(2) 

(3) 


The  non-dimensional,  scaled  variables  are  given  by 


P  = 


Re 


*  a  *  , —  * 

X  /L  ref  *  Y  =  y 


ref 


*  ,  a  ~  .  a 

u  /u  c  ,  V  =  v  /Re/u  c 
ref  ref 


A  A  A  A 

P  /P  ref»  P  =  P  /P  ref 


p  e  u  ,  L  - 
ref  ref  ref 


(4) 


(5) 


M  ref 

(See  any  standard  text  on  boundary- layer  theory,  such  as  Schlichting 
[26]  for  a  discussion  of  Reynolds  number  scaling.) 


The  usual  boundary  conditions  are 


impervious  wall 

i|>(x,o)  =  0 

(6a) 

no-slip 

u(x,o)  =  0 

(6b) 

far  field 

u(x,Y)  =  1  as  Y  +  « 

(7) 

where  the  stream 

function  is  defined  by 

3U> 

U  "  3Y 

v  =  - 
*  3x 

(8) 

The  introduction  of  the  stream  function  identically  satisfies 
continuity,  Eq.  (1). 

For  a  direct  solution  of  the  boundary- layer  equations,  an  edge 
velocity  or  pressure  coefficient  distribution  is  typically  specified. 
Inverse  solutions  of  the  boundary- layer  equations  can  be  obtained  by 
specifying  either  a  wall  shear  (or  function  of  the  wall  shear)  or  a 
displacement  thickness  distribution.  The  form  of  the  resulting  boundary 
conditions  will  be  discussed  later. 

2.2  Transformation  of  the  Boundary- Layer  Equations 

The  governing  equations  (1)  -  (3)  and  the  boundary  conditions  (6)  - 
(8)  are  transformed  using  the  modified  Levy-Lees  pseudo-self-similar 
transformation.  Once  the  equations  are  in  similarity  form,  the  FLARE 
approximation  of  Reyhner  and  Flugge-Lotz  [5]  can  be  introduced  to  handle 
thin  regions  of  separation.  A  composite  coordinate  transformation  is 
then  introduced  to  resolve  the  two  length  scales  present  in  turbulent 
boundary  layers.  The  advantages  and  disadvantages  of  these 
transformations  are  also  discussed  later. 
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2.2.1  Modified  Lew-Lees  Transformation 


The  Levy-Lees  transformation  is  well  known  and  has  proven  to  be 


very  useful  for  capturing  the  growth  of  laminar  boundary- layers .  The 


standard  Levv-Lees  transformation  (intended  for  laminar  boundary  lavers) 


is  not  suited  for  capturing  the  growth  of  turbulent  boundary  layers 


because  the  laminar  viscosity  is  so  much  smaller  than  its  turbulent 


counterpart.  For  turbulent  layers,  a  transformation  which  captures  the 


rapid  growth  of  the  turbulent  boundary  layer  and  properly  resolves  the 


wall  and  wake  regions  of  the  layer  is  needed.  Carter,  et  al.  [25]  have 


developed  such  a  transformation.  This  transformation 


(1)  captures  the  growth  of  the  turbulent  boundary  layer  by 


modifying  the  form  of  the  Levy-Lees  transformation,  and 


(2)  automatically  scales  the  inner  and  outer  regions  of  the 


boundary  layer  using  a  composite,  adaptive  coordinate 


transformation.  (This  coordinate  is  discussed  in  Section 


2.2.3) 


The  modified  Levy-Lees  variables  are  given  by 


C(x)  =  0f  pe  Ut  Ue  dX 


n  ( x ,  Y )  = 


PeUe  Y  o 

■  ±r  o'  r 

✓2  5  6 


where 


Mt  =  ue  +  (pe)ref 


For  a  laminar  boundary  layer,  £ref  =  0  and  the  modified  transf ormatic 


reverts  to  the  standard  Levy-Lees  form.  For  a  turbulent  layer,  £ref 


Emax  an<*  typically,  £ref  =  eouter’  assuming  a  constant  value  of  the  eddy 


viscosity  in  the  wake  region. 
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and 
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in  i  ,1  due  i  i 
35  ly  *  3T  '  ~)n 

e  e 


Using  the  definition  of  g,  this  yields 


In  i  _  n_  fa.  i  ■) 

35  'y  25 


(14) 


The  transformed  boundary  conditions  are  obtained  from  Eqs.  (6) 
and  (7),  yielding 

f(5,0)  =  fn(5,0)  =  0  (15a) 

fq  (5,n)  =  1  as  n  nM  (15b) 

2.2.2  FLARE 

When  solving  regions  of  separated  flow  with  the  boundary- layer 
equations  using  a  marching  technique,  stability  problems  are  encountered 
because  of  the  parabolic  nature  of  the  equations.  To  prevent  the 
solution  method  from  becoming  unstable,  the  equations  must  be  modified 
inside  the  separated  region.  A  standard  procedure  to  stabilize  the 
solution  in  thin  separated  zones  was  given  by  Reyhner  and  Flugge-Lotz 
[5]  and  is  known  as  the  FLARE  approximation.  In  the  FLARE  approach,  the 
streamwise  convection  term  in  the  boundary- layer  equations  is  neglected 
inside  the  separated  region.  This  stabilizes  the  marching  procedure, 
and  since  the  convection  terra  is  small  in  thin  regions  of  separation, 
introduces  only  minor  errors  into  the  solution.  More  complicated 
procedures,  such  as  that  devised  by  Williams  [6],  require  several 
iterations  downstream  and  upstream  through  the  separation  bubble.  While 
being  more  accurate,  these  procedures  can  only  be  used  when  the  boundary 
layer  reactaches  or  the  separated  region  reaches  an  asymptotic  state. 
Since  only  thin  regions  of  separation  (backward  flow  velocities  no  more 


* 
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than  102  of  the  edge  velocity),  are  considered,  some  of  which  do  not 
close,  the  FLARE  approach  was  chosen  here. 

To  apply  the  FLARE  approximation  to  the  untransformed  x-momentum 
equation  (2),  a  FLARE  coefficient  0  is  defined  so  that 


where 


o  3U  .  „  3u  _ e  3  ,  9uv 

0u  r—  +  V  —  =  u  - —  +  TT7  (v.  tit) 
3x  3Y  e  dx  3Y  t  3Y 


0  =  { 


1 ,  u  >  0 


0,  u  <  0 

Horton  [11],  however,  notes  that  when  using  pseudo-self -similar 
variables,  this  approach  is  too  restrictive,  since  it  affects  the  self¬ 
similar  as  well  as  the  non-similar  solutions. 

To  devise  a  less  restrictive  approximation,  the  following 
substitutions  are  defined  after  noting  the  form  of  the  convection  terms 
on  the  right-hand  side  of  Eq.  (12): 

f  =  —  =  U 

n  % 

f5  =  -v 

Substituting  these  definitions  into  Eq.  (12),  the  right-hand  side 
becomes 

*  ■ 

This  result  parallels  the  form  of  the  convection  terms  in  the 
untransformed  equation,  (2),  with  £  and  n  replacing  x  and  Y 
respectively.  Since  the  streamwise  convection  term  introduces  the 
instabilities  in  regions  of  backflow,  only  the  analogous  transformed 
term  needs  to  be  removed.  The  FLARE  variable  is  applied  only  to  the 
U  —  term,  and  Eq.  (12)  becomes 

(bf  )  +  ff  +  P(1  -  f2)  *  2?[0f  f,  -  f.f  ]  (16) 

nn  n  nn  n  n  5n  5  nn 

where 
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1,  f  >  0 

e  .  n  ~ 

u  o,  f  <  o. 

n 

Here,  when  5  is  set  to  zero  for  a  Falkner-Skan  or  starting  solution,  no 
0-terms  remain  which  would  interfere  with  obtaining  the  self-similar 
solution  that  is  sought. 

2.2.3  Adaptive  Composite  Coordinate 

The  modified  Levy-Lees  transformation  captures  the  growth  of  the 
boundary  layer  for  both  laminar  and  turbulent  flows,  but  two  length 
scales  still  exist  for  the  turbulent  case  and  the  wall  region  must  be 
adequately  resolved.  To  accomplish  this,  a  new  normal  coordinate  is 
introduced  for  the  inner  region, 

Vcf 12  ,  r - 

N.  *  - ? -  tan  (a/RC.  £  n)  (17) 

i  3  /  6  r 

V  e 

This  inner  coordinate  is  based  on  an  approximate  analytical  velocity 
profile  given  by  Whitfield  [ 27 ] . 

The  composite  coordinate  is  formed  by  matching  to  an  outer 
normal  coordinate  which  is  derived  from  the  universal  coordinate  work  of 
Clauser  [28].  The  outer  coordinate  has  the  form 

N  =  (1  -  b)  tanh  l/a  {  d[(l  -  b)  r,  /  <$*]“}  +  b  (18) 

o 

where  b  is  a  slip  velocity  at  n  =  0,  allowing  for  the  presence  of  the 
inner  layer. 

The  composite  is  given  by 
N  =  N.  +  N  -  N. , 

l  o  I/O 

where  is  the  outer  limit  of  the  inner  coordinate.  This  is  matched 

to  the  inner  limit  of  the  outer  coordinate,  giving  the  result 
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Ni/o  "  b  "  2a 


(19) 


The  composite  coordinate  maps  the  serai-infinite  line  0  <  q  <  <*>  to 


0  <  N  <  1. 

Carter  et  al.  [25]  test  the  composite  coordinate  in  this  form, 
using  a  noniterative  finite-difference  scheme,  with  satisfactory 
results.  The  coordinate  in  this  form  applied  to  the  present  scheme  has 
several  problems.  These  problems  are  due  mainly  to  the  asymptotic  form 
of  the  outer  coordinate.  As  n  approaches  infinity,  N  approaches  unity 
and  the  metric  coefficient 


This  metric  appears  in  the  denominator  of  the  coefficients 
resulting  from  the  spline  discretization  (See  section  2.7).  Even  when 
the  coordinate  N  is  taken  to  be  slightly  less  than  unity,  as  discussed 
by  Edwards  et  al.  [29],  the  division  by  small  metrics  remains  a  problem. 
To  correct  this  problem,  the  outer  coordinate,  Eq.  (18)  as  given  by 
Carter  et  al.  [25],  is  modified  so  that  n  does  not  asymptotically 
approach  infinity  as  N  approaches  unity. 

After  attempting  a  linear  and  quadratic  distribution  for  NQ,  a 
coordinate  based  on  the  inverse  hyperbolic  sine  was  adopted.  This 
coordinate,  while  maintaining  close  similarity  with  the  hyperbolic 
tangent,  does  not  approach  an  asymptote  and  so  does  not  create  problems 
with  the  metric  coefficients. 

The  inner  coordinate  is  kept  as  before, 

tan  ^  (a2h) 


>?> 

jA. 


’V  s 

-  V 


<  V 


‘.  v  \  -  '  * 
■VA 


'■.V.V.v'.n' 
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where 

^yv2  /a 
“2  ■  aVSe  Cf* £ 

(20) 

The  new 

outer  coordinate  is  chosen  to  be 

N  =  (1  -  b)  a-  sinh  ^  (a.r\)  +  b 

0  j  4 

(21) 

and  the 

following  conditions  are  enforced: 

1) 

N  =  b  at  q  =  0, 

0 

2) 

N  =  1  at  n  ■  n  » 

0  *  *00* 

3) 

t*  1.  -  0  -  b,  d  ‘/-/i* 

Condition  (3)  is  prescribed  by  forcing  the  slope  of  the  new  outer 
coordinate  to  match  the  slope  of  the  original  outer  coordinate  at  n  *  0. 

Condition  (l)  is  automatically  satisfied  by  the  form  chosen  for  NQ. 
Condition  (2)  yields  the  result 

p  a.  =  sinh  (1/a-)  (22) 

00  4  j 


and  condition  (3)  yields 


a3a4  =  (1  -  b)  dl/a/S* 


(23) 


These  equations  result  in  a  transcendental  relation  for  83,  which  is 
solved  using  the  Newton-Raphson  technique,  after  which  may  be  found 
from  Eq.  (22)  or  (23).  The  modified  composite  coordinate  is  then  given 
by 

N  =  a^  tan  ^  (a^n)  +  (1  -  b)  a^  sinh  *  (a^n)  (24) 

where  b  *  ^  a^ 


•.  v*«* 


a^>  =  determined  by  Newton -Raphson  iteration 

a  »  0.09 

d  =  0.A5 

a  =  1.25 

With  the  composite  distribution  specified,  the  corresponding  n 
distribution  can  be  determined  by  again  using  Newton-Raphson  iteration. 
The  first  and  second  order  metrics,  to  be  used  later,  are  given  by 


3N 

aj^  U  ' 

‘  b>  a3a4 

3n 

i  +  <a2n)2  y7 

+  (a4n)2 

32N 

o  3 

2ala2n 

(1  -  b)  a3a4 

an2 

[l  +  (a2n)2]2 

[l  +  (a4n)2] 

(25a) 

(25b) 


With  the  adaptive  composite  coordinate  fully  defined,  the  x- 
momentum  equation  in  (5,n)  coordinates  can  be  transformed  to  (£,N) 
coordinates.  The  derivatives  transform  according  to 


(— ) 


Nn 


(26a) 


~  *VaN^  + 


(26b) 


and 


(1 - )  a  N  (•§-)  +  N2  (- - ) 

3n2  5  nn  n  V  C 


(26c) 


The  transformed  x-momentum  equation  is  put  into  first  order  form. 


yielding 
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c 

ftl 

I 


f  -  u 

n 


U  *  T 

n 


(bx)  +  fx  +  8(1  -  u  )  ■  2£(0uu  -  xf  ). 

n  t,  e. 

Transforming  to  the  (£,  N)  plane,  Eqs.  (27)  become 

N  fM  -  u 
q  N 


N  u  =  T 
q  N 


(27a) 

(27b) 

(27c) 

(28a) 

(28b) 

(28c) 


N  (bx)N  +  fx  +  3(1  -  u)  =  25[0u(N  Uj,  +  u  )  -  x(N^fN  +  f  J I  • 

Substituting  from  Eqs.  (28a)  and  (28b),  Eq.  (28c)  is  written  as 

2  NE  Nf 

N  (bx )  +  fx  +  g(l  -  u  )  *  2E(Qu(^  x  +  u  )  -  x(^  u  +  f  )]. 

n  n  n  ^ 

ne 

The  term  Q  ^  ux  never  appears,  regardless  of  the  value  of  0,  since 

n 

it  cancels  with  the  term  of  opposite  sign  when  0  is  set  to  unity,  and 
when  the  flow  is  separated,  0*0.  Thus  the  above  equation  may  be 
simplified  to 

N  (bx)„  +  fx  +  8(1  -  u2)  =  2E(0uu_  -  xf_  +  (0  -  1)  rux]  (29) 

q  N  t,  t, 

where 

r  *  N  /N 

5  n 

To  avoid  taking  the  derivative  of  the  eddy  viscosity  ratio  b,  a  new 

variable  is  introduced  as 
* 

x  *  b  G 

where 

*  - 1 

b  =*  b  =  p  /u  . 

e 


Then,  the  final  form  of  the  governing  equations  is 


N  f  -  u 
q  n 


(30a) 
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N  u  *  b*G  (30b) 

0  n 

N  Gv,  +  b^fG  +  p(l  -  u2)  +  (1  -  0)2£rb*uG  =  £[e(u2)  -  2b  Gfr]  (30c) 

n  N  t,  £, 

Here,  three  important  items  can  be  noted: 

(1)  The  governing  equations  in  (£,N)  coordinates  can  easily  be 
reduced  to  (£,n)  coordinates  by  setting 

N  =  1 

n 

N.  -  0 

t, 

(2)  The  composite  coordinate  must  be  simplified  for  laminar  flow. 
The  inner  coordinate  and  the  slip  velocity  are  set  to  zero, 
since  these  allow  for  the  presence  of  the  turbulent  wall 
layer.  For  flows  transitioning  from  laminar  to  turbulent,  the 
inner  coordinate  should  be  "switched  on"  at  the  beginning  of 
transition  to  allow  the  coordinate  to  grow  with  the  changing 
boundary  layer. 

(3)  The  term  (1-0)  2?rb*uG  does  not  appear  unless  (fj,N) 
coordinates  are  used  and  the  flow  separates.  Since  the 
adaptive  composite  coordinate  is  only  intended  for  use  with 
attached  boundary  layers,  this  term  never  actually  contributes 
in  the  problem  considered  here.  It  is  included  for 
completeness,  however,  in  the  event  that  an  adaptive 
coordinate  is  developed  for  both  separated  and  attached  flow. 
This  would  involve  modifying  the  inner  coordinate  to 
account  for  regions  of  backflow  (see  Swafford  [30]  and 

Whitf ield  [ 31] ) . 
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2. 3  Mechul  Function  Method 

The  final  form  of  the  governing  equations  in  first-order  form  is 
given  in  the  previous  section  by  Eqs.  (30),  with  the  boundary  conditions 
fU.O)  =  u(5,0)  =  0  (31a) 

u(C.n)  =  1  as  n  -  Hoc  (31b) 

where  rioo  is  chosen  as  a  suitably  large  constant.  In  this  system,  there 
are  three  first-order  equations,  three  variables  f,  u,  and  G,  and  three 
boundary  conditions.  If  (5  is  given,  as  in  a  direct  case,  then  the 
system  of  equations  is  complete. 

For  an  inverse  solution,  however,  an  additional  boundary  condition 
is  added,  which  consists  of  specifying  wall  shear,  displacement 
thickness,  or  another  quantity.  The  streamwise  pressure  gradient 
parameter  is  unknown,  and  the  system  is  overdetermined. 

To  close  the  system,  Cebeci  and  Keller  [22]  introduce  an  additional 
differential  equation.  The  form  of  this  auxiliary  equation  is 
determined  by  letting 

0(0  3  0(5, n) 

and  requiring 


This  provides  the  additional  relation  needed  to  close  the  system  of 
equations  for  an  inverse  solution. 

While  this  formulation  of  the  mechul  function  method  was  applied  in 
Reference  22,  Kaufmaa  and  Hoffman  [23]  discovered  that,  when  applying  a 
spline  f inite-difference  scheme  (particularly  spline  S^(4,0),  [21])  to 
the  mechul  function  method,  Eq.  (32)  resulted  in  a  singular  matrix. 

This  singularity  was  due  totally  to  the  form  of  the  spline  relations. 

To  prevent  this  problem,  the  auxiliary  equation  is  changed  to 


with  Eq.  (32)  enforced  through  the  boundary  conditions. 


Transforming  Eq.  (33)  with  Eq.  (26c),  the  resulting  system  of 
equations  is  given  by 

Nn  fN  =  u 

NnuN  =  b*G 

N  G.  +  b*fG  +  p(l  -  u2)  +  (1  -  0)  r  b*uG  =  £[0(u2)r  -  2b*GfJ 
n  N  E,  E, 


(34a) 

(34b) 

(34c) 

(34d) 


with  boundary  conditions  given  by  Eqs.  (31),  (33),  and  the  specified 
inverse  condition,  to  be  discussed  in  the  following  section.  These 
relations  close  the  system  so  that  it  can  be  solved. 


2.4  Inverse  Formulations 

The  system  of  equations,  including  the  momentum  equation  in  first 
order  form  and  the  auxiliary  equation  obtained  through  the  mechul 
function  method,  is  given  by  Eqs.  (34).  Three  of  the  necessary  boundary 
conditions  are  given  by 


f(5. 

0)  =  0 

(35a) 

u(C, 

0)  =  0 

(35b) 

u(5. 

n»)  =  l 

(35c) 

with  Eq.  (32)  to  be  applied  either  at  the  wall  or  far  field,  as  needed. 
For  an  inverse  sr,1,< tion,  an  additional  equation  is  needed  which 
specifies  the  van  shear,  displacement  thickness  or  similar  quantity. 
This  section  discusses  four  possible  boundary  conditions  which  are 


considered  in  this  study.  Two  conditions  specify  a  function  of  the  wall 
shear,  while  the  remaining  two  specify  the  displacement  thickness. 

2.4.1  Transformed  Mall  Shear  Method 

The  inverse  solutions  of  Keller  and  Cebeci  [8],  Cebeci  and  Keller 
[22],  where  the  mechul  function  is  introduced,  Horton  [10],  and  Kaufman 

and  Hoffman  [23]  specify  the  self-similar  form  of  the  wall  shear,  given 

!  ( 

by  f  (where  primes  denote  differentiation  with  respect  to  q ) .  This  is 
by  far  the  simplest  of  the  inverse  boundary  conditions  when  the 
equations  are  cast  in  similarity  form.  The  boundary  condition  takes  the 
form 

C  -  s<«> 

In  first  order  form,  in  terms  of  the  variable  G,  this  becomes 

b*  Gw  =  S(0  (36) 

This  condition  with  Eqs.  (32)  and  (35)  closes  the  system  of  equations  in 
differential  equation  form.  When  the  equations  are  discretized, 
numerical  boundary  conditions  are  needed  to  complete  the  resulting  block 
tridiagonal  system.  These  additional  conditions  are  discussed  in 
Section  2.7. 

2.4.2  Skin  Friction  Coerficient  Method 

Although  Eq.  (36)  is  the  simplest  condition  to  apply,  it  is  not  the 
most  useful,  .rora  a  physical  standpoint,  specification  of  a  skin 
friction  coefficient  is  more  plausible.  Here,  a  skin  friction 
coefficient  with  one  of  two  possible  definitions  is  specified: 

(1)  Skin  friction  coefficient  based  on  the  reference  velocity  or 
the  velocity  in  the  freestream, 
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w 


1  *  ,  *x 

2  p  (uJ 


2  ' 


(2)  Skin  friction  coefficient  based  on  the  boundary- layer  edge 
velocity. 


Cf  = 
e 


2 

1  k  k 
\  <■  <-.> 


Transforming  these  definitions  yields 


2u  ,  , 

C.  Ste  =  -s  f 

f-  J2 %  ” 


(37a) 


C  =  —  f 

fe  Jll  W 


(37b) 


When  a  skin  friction  coefficient  is  specified,  however,  f  no 

■  i 

longer  remains  constant.  To  obtain  a  boundary  condition,  f  is 


w 


specified  as  a  function  of  or  ,  ug  and  £,  so  that 


f«  *  F(cf  v  5) 


where,  for  C, 


F  =  \  yi£  (Cf .  >/Re  ), 
e 


(38a) 


and  for  C, 


F  =  2  “F  (Cf 

u  ® 

e 

The  boundary  condition  takes  the  form 


(38b) 


bw  Gw  •  F  =  °* 


(39) 


*1 


This  boundary  condition  also  closes  the  system  of  equations  in 
differential  equation  form,  and  as  discussed  in  the  previous  section, 
additional  conditions  are  required  to  complete  the  system  in  discretized 


2.4.3  Hybrid  Method 

The  hybrid  formulation  was  developed  as  a  possible  method  for 
specifying  the  displacement  thickness  while  still  making  use  of  the 
shear  boundary  condition  at  the  wall.  This  hybrid  condition  is 
implemented  through  the  integral  momentum  equation,  in  similarity  form, 
as  a  wall  boundary  condition. 

The  x-momentum  equation  (12),  with  b  set  to  unity  for  laminar  flow 
only,  is  integrated  across  the  boundary  layer,  yielding 


f”  =  (3  +  1)  Q  +  P<5*  +  2  5  || 


where 


o  ■  q/  f  (l  -  f  )  dn 


6  *  0 /  (l  -  f  )  dq 


H  =  5  /Q 


(41a) 


(41b) 


(41c) 


The  similarity  form  of  the  momentum  thickness,  0  is  eliminated  from 
Eq.  (40)  in  favor  of  the  shape  factor  H,  which  tends  to  vary  at  a  slower 
rate.  The  resulting  equation  is  given  by 


fw'  =  5*  (Sr^ +  3)  +  2S  d^  (5*/H) 


Eq.  (42)  permits  the  specification  of  5  while  maintaining  the 


(42) 


fjj 

#J 


inverse  shear  condition  at  the  wall.  This  condition  introduces  0 


dependence  as  well  as  the  shape  factor,  H,  which  varies  with  £.  To 


handle  this  problem,  the  shape  factor  is  assumed  to  be  a  function  of  3 


alone.  This  is  strictly  true  for  self-similar  solutions  only,  but  is  a 


better  approximation  than  simply  neglecting  the  dependence  on  the  shape 


factor  in  Eq.  (42), 


2.4.4  Displacement  Thickness  Method 


A  more  direct  method  of  specifying  the  displacement  thickness 


contrast  to  the  hybrid  method,  involves  obtaining  a  boundary 


condition  containing  <S  to  be  enforced  at  the  far  field.  The 


standard  definition  of  displacement  thickness  is  given  as 


«  co  u  . 

6  -  (  1  -  ->  dy. 


Stretching  the  normal  coordinate  and  transforming  to  similarity 


variables,  Eq.  (43)  becomes 


*  —  /2E  ''<*>  ,  ' 

6  /Re  =  — *  q/  ( 1  -  f  )  dn 

e 


The  transformed  displacement  thickness  has  been  previously  defined 


to  be 


6  =  0/  (1  -  f  )  dn 


Integrating  across  the  boundary  layer,  Eq.  (41b)  becomes 


6  =  n  ~  f 


as  n  -*■  n 


f  =  n  *  $ 


If  the  untransformed  displacement  thickness  is  specified,  the 


substitution  from  Eq.  (44)  yields 


AlJil  -J*  ,  -J*  -T  .  ^  .A  .A.,  wfV-  -*  J 


This  condition,  as  opposed  to  the  conditions  obtained  in  the  previous 
three  subsections,  is  applied  at  the  far  field,  where  n  =  ria>. 


2.5  Direct  Formulation 

The  system  of  equations  for  obtaining  an  inverse  solution  of  the 
boundary- layer  equations,  as  discussed  in  this  study,  Eqs.  (34),  can 
also  be  used  to  solve  for  a  direct:  boundary- layer  solution.  For  a 
standard  direct  boundary- layer  solution,  the  streamwise  pressure 
gradient  is  given,  either  directly  or  indirectly  by  specifying  ue  or  Cp 
and  p  is  treated  as  a  fixed  parameter  at  each  streamwise  station. 

For  the  system  of  equations  considered  here,  p  is  considered  as  a 
variable  to  be  obtained  during  the  solution  process.  Thus,  to  solve  a 
direct  case  using  this  system,  an  additional  boundary  condition  must  be 
specified  which  sets  p  at  each  streamwise  station.  The  boundary 
conditions  to  be  set  are  then 


f(5. 

0)  =  0 

(47a) 

u(5. 

0)  =  0 

(47b) 

u(5. 

O  =  i 

(47c) 

3(5. 

O  =  3(5) 

(47d) 

Here,  Eq.  (47d)  takes  the  place  of  the  inverse  boundary  condition 
discussed  in  the  previous  section. 

A  direct  solution  is  typically  obtained  by  specifying  the  edge 
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Given  Cp,  the  normalized  edge  velocity  is  easily  determined  from  the 

above  relation.  To  obtain  a  8  distribution,  Lagrange  polynomials  are 

used  to  form  backward  differences  (See  Appendix  A).  For  laminar  flows, 

the  £  and  8  distributions  can  be  determined  prior  to  beginning  the 

solution  iteration  process.  For  turbulent  flows,  however,  £  and  8  must 

be  updated  after  each  iteration  because  of  the  presence  of  p  in  the 

e 

integral  for  £.  Thus,  these  quantities  can  only  be  computed  as  the 
solution  proceeds  downstream. 

2.6  Eddy  Viscosity  Model 

An  algebraic,  two-piece  eddy  viscosity  is  chosen  as  the  turbulence 
model  for  the  present  study.  This  type  of  model  is  based  on  the  original 
assumption  of  Boussinesq  that  the  apparent  turbulent  shearing  stresses 
might  be  related  to  the  rate  of  mean  strain  through  a  scalar  turbulent 
viscosity.  Algebraic  models,  while  simple  in  form  and  application,  have 
been  shown  to  yield  satisfactory  results  for  most  boundary- layer 
problems  [32].  The  eddy  viscosity  is  also  compatible  with  the  modified 
Levy-Lees  transformation.  The  eddy  viscosity  is  calculated  in  kinematic 
form,  so  that  the  effective  turbulent  viscosity  is  given  by 

“t  =  “m  +  C£ 

where  p  =  dimensionless  molecular  viscosity 
m 

e  =  dimensionless  kinematic  eddy  viscosity 
The  eddy  viscosity  model  used  here  was  developed  by  Baldwin  and 
Lomax  [24].  This  model,  while  using  the  standard  Prandtl-van  Driest 
formulation  of  Cebeci  and  Smith  [33]  for  the  inner  region,  replaces  the 
usual  Clauser  formulation  with  an  outer  equation  which  does  not  require 
the  determination  of  the  displacement  thickness.  Baldwin  and  Lomax 
maintain  that  this  is  advantageous  for  computation  of  separated  flows. 


-j 


Aka  tla 


where 


wAKE  =  r  Y  ,, 2  /v  the  smaller.  (50a) 

WA“’  lwk  yhax  ^)IF/FMAX 

In  Eqs.  (50),  is  obtained  from 

F(y)  =  y|u)|  [1  -  exp(-y+/A+)]  (51a) 

where  y,^  is  the  y  value  at  the  maximum  F(y).  The  velocity  difference 
Upjp  is  defined  as 

“DIF  =  (  Vu2  +  v2  +  w2)^  -  (  Vu2  +  v2  +  w2)min  . 

In  two  dimensions,  with  the  second  term  taken  as  zero,  this  reduces  to 
Udif  =  (  vu  +  v  )MAX  .  (51b) 

Also  in  Eq.  (49),  the  Klebanoff  intermittency  factor  is  given  by 


F^  (y)  =  (1  +  5.5  (SmV)  '! 
KLEB  yMAX 


(51c) 


The  required  constants  are 


0.0168 


C^p  —  1.6 

CKLEB  =  °-3 

C^k  —  0.25. 

For  simplicity  in  this  study,  the  Klebanoff  intermittency  factor  is 
taken  to  be  unity  across  the  layer.  The  above  equations  then  yield  an 
algebraic  eddy  viscosity  with  a  constant  value,  eQ,  in  the  outer  region. 
This  value  is  taken  as  eref. 

To  apply  this  turbulence  model  to  the  present  problem,  Eqs.  (48)- 


(51)  must  be  nondimensionalized  and  transformed  to  the  modified  Levy- 
Lees  variables.  See  Appendix  B  for  a  complete  list  of  all  transformed 
quantities  related  to  the  turbulence  model. 
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2.7  Discretization 

With  the  analytical  form  of  the  governing  equations,  boundary 
conditions,  and  turbulence  model  fully  defined,  the  discrete  form  of  the 
equations  is  next  considered.  The  equations  are  discretized  using 
fourth-order  accurate  splines  to  approximate  the  normal  derivatives  and 
backward  differences  to  approximate  streamwise  derivatives.  This  yields 
a  fully  implicit  scheme,  which  is  desirable  for  stability  and  accuracy. 


2.7.1  Backward  Differences 

The  system  of  equations,  Eqs.  (34),  is  first  written  in  index  form 
at  a  nodal  point  as 


N  .  .  f . .  =  u. . 

HiJ  ij  ij 

'  * 

N  . .  u. .  =  b. .  G. . 

HJ  ij  ij  ij 


(52a) 

(52b) 


N  . .  g! .  +  b*.  f. .  G. .  +  0. .(1  -  u2.)  +  (1 
mj  ij  ij  ij  ij  ij  ij 

=  e.  [0. .  (u2.)r  -  2  b*.  G. .  f  .] 
i  ij  ij  5  ij  ij  5iJ 


0..)  r..  b..  u..  G.. 
ij  ij  ij  ij  iJ 


(52c) 


p.j  =  0  (52d) 

where  primes  refer  to  derivatives  with  respect  to  N,  and  i  and  j  are  the 
indices  referring  to  streamwise  and  normal  location  of  the  node, 
respectively,  in  the  computational  domain.  The  final  discretized  system 
of  equations  will  be  solved  over  a  domain  where  1  <  i  <  M+l  and 
1  <  j  <  K+l. 

The  ^-derivatives  occurring  in  Eq.  (52c)  are  discretized  using 
variable  stepsize  backward  (upwind)  differences  derived  from  Lagrange 
polynomials.  In  general,  the  backward  difference  expression  has  the 
form 


/V' 


where  a=-(b+c) 


and 


b 

ii 

> 

>— > 

/Bi 

c 

I 

=  a2 

/B2 

1) 

for  ; 

3-point 

backward  differences: 

f 

A1 

"  5i 

^i-2  ’ 

B1  " 

ut-i  -  V 

(5i-l  "  5i-2)’ 

A2 

=  5i 

-  *l-v 

B2  * 

(5i-2  -  5i> 

^i-2  *  5i-l^’ 

2) 

for  ; 

3-point 

backward  differences: 

» 

A1 

-  ; 

L, 

B1 

"  ^i-l  • 

«i  ’ 

c 

*  0. 

Applying  Eq.  (53)  to  the  ^-derivatives ,  Eq.  (52c)  becomes 


'  *  7  _  * 

i  ..G..  +  b . .  f .  .  G..  +  p. .  (1  -  u.T)  +  (1  -  9 . . )r . .  b. .  u.  .  G 
mj  1J  ij  ij  ij  pij  ij'  ij'  ij  ij  ij  ij 


(54) 


=  C.[9  (au?  +  bu2  +  cu2  ).  -  2b*.  G.  .  (af.  +  bf.  .  +  cf.  ,).]. 
x  ij  x  i-l  i-2  j  xj  ij  i  x-1  i-2 


With  all  5 -derivatives  eliminated  in  the  momentum  equation,  the  system 
is  ready  to  be  discretized  in  the  normal  direction. 


2.7.2  Splines 

The  normal  derivatives  in  Eqs.  (52a),  (52b),  (52d)  and  (54)  are 
approximated  using  two  fourth  order  spline  formulas,  S^(4,0)  and 
S^(4,0),  derived  by  Rubin  and  Khosla  [21],  To  begin,  the  followirg 
spline  derivative  approximations  are  introduced: 


■? 


a 


3 

S 


1 


ft 


S 


SI 


1 


a 


^  =  @x 


lp  =  gx 


Then  Eqs.  (52a),  (52b),  (54)  and  (52d)  are  written  as 


n!  .  lf .  =  u. . 
ij  ij  ij 


(56a) 


N.  .  .  =  b.  .  G.  . 

ij  iJ  ij  iJ 


(56b) 


n!  .  «.G.  «  p..(u2.  -  1)  -  b? .  f..  G..  +  (0..  -  l)r..  b*.  u..  G.  . 
ij  ij  ij  ij  ij  ij  ij  ij  ij  ij  ij  ij 


(56c) 


+  5i  [0lj  <a“i  +  bVl  +  “1-2(1  -  2  b!j  G1J  (a£i  +  bfi-l  +  cfi-2)j 1 


L7 .  -  0 

ij 


(  56d) 


where  primes  now  denote  differentiation  with  respect  to  n*  Collecting 


terms,  Eq.  (56c)  becomes 


n! .  1°.  *  0. .  C.  u2.  +  C.  b*.  f . .  G. .  +  3. .  (u2.  -  1) 

ij  ij  ij  ij  2  xj  ij  ij  ij  ij 


+  (0. .  -  1 )r .  .  b*.  u. .  G. .  +  C.  b* .  G. .  +  0. .  C. 
ij  iJ  IJ  ij  ij  3_  ij  xj  ij  4_ 


where 


Cl..  “  ^i 
ij 


C2  *  -  (2aEi  +  1) 

ij 


Spline  S  (4,0),  defined  by. 


2.g  +  (1  +  a)2  l?  +  a2  2.8  , 

J  +  l  J  J-l 


(58) 


2  l  rl+2a  .  a  -  1  , ,  .  s3  2,„  .  v 

TT^  ”  [” —  gj+i  +  ~T~  (1  +  o)  gj  '  °  (2  +  a)gj-i] 


with 


h  .  =>  AN  . 
J  J 


and 


o  *  a .  *  h  .  /h . 
J  J  +  l  J 


is  introduced  to  eliminate  the  spline  derivative  approximations  2.^,  £u, 
and  2.^  from  Eqs.  (56a),  (56b)  and  (57).  Considering  each  of  these 
equations  in  turn  yields  (with  i  subscript  understood) 


<’2(7)j-l  +  3l(7)j  +  ^j+l  -  °3  fj-l  +  °4fj  +  °5  £j+l'  (59a) 


*  *  * 

/bG\  ,  /b  Ga  ,  b  Gx  rcru.' 

°2  +  °1  +  (7_)J+1  *  °3uj -1  +  °4uj  +  °5  uj  +  rl55b' 


and,  defining  the  right-hand-side  of  Eq.  (57)  as  RHS^j. 


n!  .  lG  .  =  RHS .  .  , 
ij  ij  ij 


the  final  first-order  spline  relation  is, 


a_  (— j— )  .  .  +  o.(— — ) .  +  (— i — ).,,  *  a.  G.  +  a.  G.  +  o.  G  ... 
2  m  J’1  1  M  J  XI  J  +  1  3  j-l  4  j  5  j+l 


(59c: 


See  Appendix  C  for  a  list  of  all  spline  coefficients. 


Spline  S  (4,0),  defined  by, 


I  3  2  2 

o  •♦•o-l  .  g  a  +  4q  +  4o  -f  1  g  1  +  q  -  q  g 

12q  j+1  12o  j  12  Lj-1 


1  rl  1  +  q  A  . 

.2  q  8j+l  '  q  8j  8j-l'] 


is  introduced  to  eliminate  the  spline  derivative  approximation  from 
Eq.  (56d),  yielding 

Sj-1  +  °7  8j  +  °6  Bj+1  '  <6l> 

Eqs.  (59)  and  (60)  form  a  nonlinear  4X4  block  tridiagonal  system  for 
each  streamwise  station,  i.  The  boundary  conditions  to  complete  the 
system  are  discussed  in  the  following  section. 


2.8  Boundary  Conditions 

The  system  of  equations  obtained  above,  Eqs.  (59)  and  (60),  apply 
to  the  interior  points  of  the  4X4  block  tridiagonal  system.  To  close 
the  system,  four  boundary  conditions  must  be  supplied  for  both  the  wall 
and  far  field.  When  an  insufficient  number  of  physical  conditions  is 
provided,  numerical  conditions  obtained  from  two-point  spline  relations 
are  added  to  close  the  system.  In  the  following  subsections,  the 
boundary  conditions  for  the  shear/hybrid  methods,  displacement  thickness 
method  and  direct  method  are  discussed. 


2.8.1  Shear  and  Hvbrid  Method  Conditions 

The  following  Dirichiet  boundary  conditions  have  been  specified 


from  physical  constraints  (with  i  subscript  understood): 


(47a) 


(47b) 


(47c) 


For  the  transformed  wall  shear  case,  the  condition  given  by  Eq.  (36) 


E 


is  also  specified.  When  a  (C,  /Re)  distribution  is  specified,  Eq.  (36) 


is  replaced  by  Eq.  (39),  so  that  (with  i  subscripts  understood), 


b*  G.  -  ^  /25  (Cf  /Re)  =  0  . 

e 


To  specify  (C^  /Re),  Eq.  (39)  is  again  used,  but  with  F  given  by 


Eq.  (38b),  yielding 


b*  Gi  '  \  2  (Cf  *  0  * 

U  oo 

e 


Finally,  for  the  hybrid  method,  Eq.  (36)  is  replaced  by  Eq.  (42),  so 


*  -*  Pi  +  1  j  i* 

bi  Gi  -  4  (-T-  +  h’  +  24  3?  (h)  *  °- 


For  any  of  the  shear  conditions  together  with  the  derivative 


condition  on  8,  Eq.  (32), 


PN  U, N)  =  0 


specified,  the  system  is  complete  at  the  wall.  To  implement  Eq.  (32),  a 


two-point  spline  relation,  given  by  Rubin  and  Khosla  [21],  is  used. 


The  fourth-order  accurate  two-point  spline  relation  to  be  used 


throughout  this  investigation  is  given  by 


h-  h, 

e  -  e  -  — -  ( 2,®  +  2.®)  +  —  (I®  -  Lg)  =  0 

g2  gl  2  U2  V  12  U2  V  U' 
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When  g  =  8,  Eq.  (62)  becomes 


^2  1  '  ^2  11  11 
s2  -  Pl  -  r  <s2  +  +  ti  <-h  -  h  >  ■  °- 

Upon  substit  :ion  of  Eqs.  (32)  and  (34d),  the  preceding  relation  reduces 


P2  -  PL  -  0.  (63) 

With  regard  to  the  far  field  conditions,  the  spline  relation 
Eq.  (62)  must  be  used  three  times  since  only  one  physical  constraint, 

Eq.  (47c)  exists.  The  numerical  conditions  used  to  close  the  system  are 


fK+l  fK  ‘  2  ^fK+l  +  V  +  12  <fK+l  fK  J  ”  0  <64a) 


“KH  -  \  ‘  (UK+1  +  UK)  +  ^  (uK+l  '  =  0  <64b) 


hj/ii  i  i  hy .  i  ' '  '  t 

GK+1  "  Gk  2~  (GK+1  +  V  +  12  (GK+1  "  GK  ^  =  °  ^64c^ 

where  primes  denote  differentiation  with  respect  to  N.  Since  G  and  its 
normal  derivatives  approach  zero  exponentially  rapidly  near  the 

it  it 

boundary- layer  edge,  the  difference  (G^+^  -  Gj,  )  in  Eq.  (64c)  is 
expected  to  be  small  and  hence  is  neglected,  leaving 


"v+l  i  1 

GK+1  "  GK  2~  (GK+1  +  GK)  =  °* 


(64d) 


The  first-order  derivatives  can  be  eliminated  by  substituting  from 
Eqs.  (28),  To  obtain  the  second  derivative  of  f  with  respect  to  N,  note 
that,  from  Eqs.  (27), 


f  =  t  =  b  G 

nn 


Transforming  to  (£,N)  coordinates  with  Eq.  (26c)  and  substituting  from 
Eq.  (27a)  yields 


f  I 

f 


* 

(b  G 


nn  \ 
-  — u) 


(65) 


Similarly,  for  the  second  derivative  of  u  with  respect  to  N,  we 
begin  with  the  following  identity  in  (C,  n)  variables: 


u  h  t  =  ( b  G ) 

nn  n  n 


Normally,  this  relation  would  have  to  be  expanded  using  the  product 
rule,  but  since  Eq.  (64b)  is  applied  at  the  far  field,  b  equals  unity 
because  the  eddy  viscosity  is  constant  in  the  wake  region  of  the 
turbulent  boundary  layer.  The  relation  for  u  then  becomes 


N 


II  1  f  L1 

U  =  "T  (NnG  '  G) 
n2  n  Nn 
n 


Here,  and  in  Eq.  (64d),  the  momentum  equation  is  used  to  eliminate  G' . 
Upon  making  all  necessary  substitutions,  Eqs.  (64)  become 


K+l 


.  .  Vl  ,S<+1 

PK+1  "K+l  +  PR'hc  +  12  2  M'  2 

nk+i  nk 


=  0 


(66a 


\+l 

“K+l  "  \  "  2 


N, 


K+l 


GK  ^K+l  1 

+  4)  +  44  [-T7  (RHS 

N  N 

WK  K+l 


i  i 


K+l 


-  (RHS  -  4-  G)v3  =  o 

V  N 


(66b 


tiS 


a 


where. 


=  i_  .Vi  ^k_  ^c+k 
Pk  “  N  12  n’2  "  2 

K  K 


_1 _  ,*^+1  nk-h  Nc+l, 

pR+1  =  n'  12  n'2  2 

K+l  K+l 


Summarizing,  the  wall  boundary  conditions  used  are 


fi  =  o 


“i-o 


P2  -  *1  "  0 


and,  for  specified, 


bL  G:  -  S  =  0 


or,  for  specified, 
e 


b*  Gl  -  \  /25  (Cf  /Re) 
e 

or,  for  specified. 


bt  G1  •  I  ^  <C£  ^  ■ 

u  00 


The  far  field  boundary  conditions  are 
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,  .  Nc+l  ,GK+1 

PK+1  'hc+l  +  PK\  12  ^  1 2 

N„  .  , 


(67h ) 


1+1  ,GK+1 


-1  r  1 


^  ’  X  /  fw  '  X  ,  Cv.  \  .  IV  r  1,  r  1  y  _  ..p  IV  V 

\+l  "  \  '  2  (  12  '2  (RHS  ‘  '  G)K+1 


nk+i  nk 


-4t  (RHS  -  — -  G),J  =  0 
N  ~  N  ' 


(671) 


G,,  -  G„ 

N+l  N 


\+1  RHS  RHS 
^  1  (— r^1  +  —A)  =  o 


( 6  7  j  ) 


These  are  the  conditions  used  when  a  form  of  the  wall  shear  is 
specified,  or  the  hybrid  method  is  used  with  the  displacement  thickness 
specified.  For  the  displacement  thickness  method  or  a  direct  solution, 
other  modifications  must  be  made. 


2.8.2  Displacement  Thickness  Method  Conditions 

When  the  displacement  thickness  method  is  used,  the  boundary 
conditions  described  in  the  previous  subsection  must  be  altered,  both  at 
the  wall  and  far  field.  The  inverse  condition  for  the  displacement 
thickness  method  is  implemented  at  the  boundary-layer  edge.  Thus,  Eq. 
(36)  must  be  replaced  by  a  two-point  spline  relation  and  one  spline 
condition  at  the  wall  must  be  discarded  to  allow  for  the  condition  on 


The  two-point  spline  relation  chosen  for  the  wall  condition  is 


^9  !  I  ^9  r  r  ii 

f2  '  fi  -  -f  (f2  +  +  it  (f2  '  fi  )  =  0 


This  condition  is  identical  to  Eq.  (64a)  applied  at  the  far  field  for 


the  shear  case.  After  making  the  appropriate  substitutions,  the 
conditions  are  given  by 


fK+1  -  nK+1  ♦  *  0  W 

depending  upon  whether  the  transformed  or  untransformed  displacement 
thickness  is  specified. 

2.8.3  Direct  Method  Boundary  Conditions 

For  the  direct  case,  a  slightly  different  set  of  boundary 
conditions  is  used.  At  the  wall,  the  conditions  are  identical  to  those 
used  for  the  displacement  thickness  method.  These  conditions  are  Eqs. 
(47a),  (47b),  (63),  t.id  (69). 

At  the  boundary- layer  edge,  two  physical  conditions  are  given, 

Eq.  (47c)  which  is  always  specified  regardless  of  method,  and  Eq.  (47d), 
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'  * 

« 

■? 


Two-point  spline  relations  are  required  to  close  the  system.  The 
obvious  choices  are  Eqs.  (66a)  and  (66b),  since  these  complete  the 
system  and  do  not  neglect  any  terms,  as  was  done  in  Eq.  (66c). 

It  is  apparent,  considering  the  boundary  conditions  discussed  in 
the  previous  subsections,  that  only  minor  modifications  need  to  be  made 
to  the  boundary  conditions  even  though  a  range  of  methods  is  treated 
here.  This  is  quite  surprising  and  also  very  helpful,  since  it  allows 
the  various  methods  to  be  combined  into  a  single  computer  code  with 
little  difficulty. 

It  should  also  be  noted,  however,  that  while  only  minor  changes  to 
the  boundary  conditions  are  required,  very  few  other  possibilities  exist 
when  using  the  two-point  spline  relations.  The  conditions  obtained  in 
the  above  sections  consist  of  all  available  conditions  for  the  current 
problem  without  further  differentiating  the  momentum  equation,  which 
would  greatly  increase  the  complexity  of  the  boundary  conditions. 


2.9  Linearization 


The  system  of  block  tridiagonal  equations  and  boundary  conditions 
discussed  in  the  above  sections  is  nonlinear.  The  system  is  linearized 
and  solved  using  Newton’s  method,  where  a  variable  at  the  most  recent 
iteration,  (n+l),  is  equal  to  the  variable  at  the  last  iteration,  (n), 
plus  a  small  correction.  Then,  the  variables  at  the  (n+l)  iteration  can 
be  written  as 

f(n+D  _  f(n)  +  ^(n) 


(n+L)  _  (n)  .  (n) 

u  =  u  +  6u 


G(n+1)  _  Q(n)  +  5G(n) 


3(n+l)  ,  g(n)  +  6jJ(n) 


(70) 


r  V  y-w- y  TVT7 


IJU  m  I  j  mwu 


y 
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These  relations  are  substituted  into  the  nonlinear  equations  at  the 
(n+1)  iteration  and  all  terms  which  are  quadratic  in  the  corrections  are 
neglected.  This  yields  a  linear  system  of  equations  which  can  be  solved 
for  the  corrections  5f,  5u,  <SG,  and  63  using  a  lower-upper  (L-U) 
decomposition  method.  The  following  subsections  describe  the 
linearization  of  the  eddy  viscosity,  the  interior  spline  relations  and 
the  boundary  conditions. 

2.9.1  Eddv  Viscosity 

Before  beginning  the  linearization  of  the  equations,  we  considered 
the  turbulence  model  formulas.  In  simplified  form,  the  eddy-viscosity 
relations  can  be  written  as 

£i  *  C1  ^  I f  I  *  (71a) 

e  =  C_.  (71b) 

o  Z 

Here,  the  inner  equation  is  a  function  of  n  and  f *  (primes  denoting  r|" 
derivatives),  while  the  outer  equation  is  a  constant  when  the  Klebanoff 
intermittency  factor  is  set  to  unity.  Thus,  to  maintain  consistency  in 
the  solution,  it  is  obvious  that  the  eddy-viscosity  equation  must  also 
be  linearized. 

Upon  substituting  for  f'  ,  Eq.  (71a)  becomes 

c.  =  C^n)  b*  ! G |  ,  (72) 

where  b  can  be  taken  outside  the  absolute  value  sign  since  its  value  is 

£ 

always  positive.  Since  b  is  also  a  function  of  e^,  the  value  of  is 
determined  by  solving  a  quadratic  equation,  yielding 

e.  -  \  [(1  +  4C1(n)  (1  +  eQ)  |G| ) 1/2  -  1]  . 

To  eliminate  the  absolute  value  of  G,  Eq.  (72)  is  recast  as 
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*  (-1)S  C1(n)b  G, 


where  s  »  1  for  G  <  0,  and 


s  =  2  for  G  >  0. 


Writing  in  Newton  iterate  form, 

(n+1)  _  (n)  (n) 

e .  =  e .  +  <5e .  , 


Eq.  (73)  may  be  substituted  at  (n)  and  (n+1)  iterations  and  the  result 
solved  for  <Sefn^  . 

l 

Se!n)  =  (-i)s  cl  (n)  (b*G)(n+1)  -  (-l)s  C^n)  (b*G)(n)  . 

Upon  substitution  for  the  (n+1)  Newton  iterates,  as  in  Eq.  (70),  the 
above  equation  may  be  solved  for  the  correction  of  the  solution  variable 
G,  yielding 


,.<■»  .  <;*>  5G(n, 

1  (2e.^  +  1) 

i 


(75a) 


or,  in  a  more  convenient  form. 


(-i)s  c  (n)  b*  (n)  ,  , 

- i - 7-T -  6G(n  ■ 

e,  (n) 


(75b) 


l  +  <r+*r> 


2.9.2  Interior  Equations 

The  interior  block  tridiagonal  equations  are  given  by  Eqs.  (59)  and 
(61).  Two  of  these  equations  contain  nonlinear  terms,  Eqs.  (59b)  and 
(59c).  Each  of  the  nonlinear  terms  is  linearized  separately,  then 
substituted  into  the  appropriate  equations.  The  remaining  equations 
are  then  linearized  and  put  into  correction  form. 

*  * 

A  recurring  term  in  both  nonlinear  equations  is  b^G^.  First,  b j , 
is  linearized  using  the  binomial  expansion  theorem,  yielding 


For  the  inner  region,  6e  is  given  by  Eq.  (75)  and  for  the  outer  region, 
6e  us  equal  to  zero.  Upon  linearization  and  substitution  of  Eqs.  (73), 

•ft 

(75b)  and  (76),  the  term  (b^G^)  becomes  (with  n  iterate  superscripts 
understood) , 

(b*Gj(n+l)  =  b*[Gj  +  (1  -  e j )  5G  ]  (77) 


■z - ; — r  for  inner  region 

2e.  +  1  ° 

l 

e .  = 

^  0  for  outer  region  . 

The  remaining  terms  in  linearized  form  are: 

(b*f .G.)(n+1)  =  b!f.  [G.  +  (1  -  e.)  6G . ]  +  b*  G.  6f.  (78a) 

JJJ  JJJ  J  J  JJJ 

(b*  u.  G.)(n+1)  =  b*u.  [G.  +  (1  -  e.)  6G.]  +  b*G.6u.  (78b) 

JJJ  JJJ  JJ  JJJ 


,  2x(n+l)  2  „  . 

(u. )  =  u.  +  2u.6u. 

J  JJJ 


(78c) 


9  (n+1)  _  „ 

[8.(u?  -  1)]  =  g.(u.  -  1)  +■  28. u.  6u.  +  (u.  -  1)  68.  (78d) 

J  J  J  J  J  J  J  J  J 

Upon  substitution  of  these  terms  into  Eq.  (59c),  the  equation  is 

solved  for  all  correction  terms,  with  all  known  quantities  moved  to  the 

right-hand-side.  The  final,  linear  tridiagonal  result  is  given  by 


[A6f  +  B6u  +  D5G  +  G6S ] . 


-7  [A6f  +  B5u  +  E5G  +  G5@ ] . 


(A6f  +  B6u  +  F6G  +  G6Bh  +  1 
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°2  N  °3 

7-  iTi  *  (t2  -  r  r> G  +  ec4>j-i 

Nj-i  j  1 


a.  o, 

4  IT.  MT,  -  J-  7)  G  +  GC.J  . 

N.  1  2  Vl  4  J 

J 


4CTl+(T2'^°5>CtSC4'j+I 

Nj+1  J 


See  Appendix  C  for  a  list  of  all  block  tridiagonal  coefficients. 

The  remaining  three  tridiagonal  spline  relations  in  linearized  form 


O-  C -  o.  o.  0»  . 

6f  -  -4  fiu).  .  +  (4  fif  -  4  fiu).  +  (4  6f  -  4  6u).+1 

hj  N  J'1  hj  N  J  hj  N  J+1 

=  (4  u  -  4  f).  ,  +  (4  u  -  ~  f>.  +  (4  -  ?r  f)i+i 

N  hj  J_1  N  hj  J  N  hj  J+1* 


(81a) 


(7  5u  4  b  (1  -  e)  6G)  +  (4  -  4  b  (1  -  e)6G). 

j  N  J'1  j  N  J 


+  au  -  |r  a  -  so)  -  (lib’s  4  u) 

J  J  J 

°1  *  °A  v*  0, 

+  (7  b  G  ‘  “  u)j +  (7 G  ‘  “  u,j+i  • 

N  j  J  N  j  J 


(81b) 


6Bj-l  +  °7  5Bj  +  a6  6Bj+l  =  '  (3j-l  +  °7  Bj  +  °6  3j+l } *  (81c) 

(See  Appendix  C  for  a  summary  of  all  spline  coefficients.) 


2.9.3  Boundary  Conditions 

Since  in  general,  the  boundary  conditions  are  linearized  by  the 
same  procedure  discussed  in  the  previous  subsections,  a  detailed 


linearization  of  the  boundary  conditions  is  not  provided  here.  Several 


special  cases  do  exist,  however,  for  the  various  "inverse"  conditions. 
These  special  cases  are  examined  here,  followed  by  a  brief  summary  of 
all  boundary  conditions  in  linearized  correction  form. 

(1)  Skin  friction  and  displacement  thickness  methods 

For  the  transformed  wall  shear  case,  the  correction  form  of 
Eq .  ( 36  )  is 

A 

(1  -  ex)  b“  SGj  -  0  , 

since  S(5)  is  a  given  constant  at  each  streamwise  station.  When  either 

or  is  specified,  however,  the  function  F  in  Eq.  (39)  varies  as 
6  °° 

the  solution  converges,  since  it  is  dependent  upon  £  and  u^  (for 

oo 

case  only).  Thus,  F  must  be  put  into  linearized  form.  Upon 
linearization,  Eq.  (39)  becomes 


(1  -  e.)  b,  6G.  -  (F_  6C,  +  Fu  6u  +  Fr  6.)  *  0,  (82) 

111  f  ee  5  5 

where  we  set  to  zero  since  is  specified.  The  partial  derivatives 

of  F  are  easily  found  to  be: 

(a)  for  Cc 


F  =  0  , 
u 

e 


(83a) 


F.  -  \  ~ZZ  (C,  /lie)  , 
5  2  /25  e 


(83b) 


b)  and  for  C, 


c  , — 

~  (Cf  /Re), 

u  00 


(83c) 


F  =  — 1 -  (I  -  p)  (c  /rI), 


■5  r-  2  2 

m  ue 


(83d) 


Returning  to  Eq.  (82)  and  substituting  the  previous  results  for  6£ 
and  6ue,  the  boundary  condition  takes  the  general  form, 

(1  -  e,)  b*  6G,  -  X.  68,  =  0  (Rl) 


1  (1  -  tx  t2)  £  4(1  -  tL  t2)  ue 

For  the  displacement  thickness  method,  a  similar  procedure  is 
followed,  except  that  the  resulting  boundary  condition  contains  6f  and 
68  insteady  of  6G  and  68.  The  final  form  is  given  by 

5fK+l  +'X1  68K+1  =  nK+l  "  fK+l  ^  ^88a 

where 


F  -  (f|)  +  u  (2£)"3/2  O  -  1)  (6*^), 
£  9£  e 


(2C)~1/2  (6*/Re), 


the  ^'derivative  j.s  found  by  backward  differences,  and  X^  is  given 
above.  If  6  is  specified,  Eq.  (88a)  is  replaced  by 


5fK+l  ~  nK+l  '  fK+l  "  5 


(2)  Hybrid  method 

Upon  discretizing  the  ^-derivative,  Eq.  (47)  becomes  at  j  =  1 

~ *  “ a 

bX  •  K  +  V  -  ^b  +  b  ■ 0  <9o: 

i  i  i-l  i-2 

Here,  the  shape  factor  must  be  linearized  along  with  the  solution 


i 


variables.  To  accomplish  this,  H  is  assumed  to  be  a  function  of  3 


alone,  so  that 


H(n+1)  =  H(n)  +  5H(n)  _  R(n)  +  6p(n) 


Strictly,  this  is  true  only  for  self-similar  flows,  but  the 


approximation  is  used  here  to  include,  at  least  partially,  the 


dependence  of  the  boundary  condition  on  the  variation  of  H. 


Introducing  the  Newton  iterates  into  Eq.  (90)  and  linearizing 


yields,  in  correction  form, 


(1  -  ex)  b1  6  G1  +  X2  6p  =  P 


(91a) 


where,  with  j=l  understood. 


-*  Hi  +  1  I  dH 

X2  *  -  5i  l-V-  -  3  <Si  +  25iS  +  l>  Oi1 

i  H. 

i 


(91b) 


'*  +  25  a  +  1 

Pl'5i  - 

l 


“*  ** 

6 .  .  6 .  _  j  A 

+  3.  +  25.  b(~— •)  +  cfcj^)  -  b.G.  (91c) 

1  1  Hi-1  Hi-2  1  1 


The  derivative  (^r)^  is  obtained  by  Lagrange  backward  differences  past 


the  second  streamwise  station.  For  1  <  i  <  2,  the  derivative  is  given 


its  self -similar  value  based  on  the  current  value  of  3* 


(3)  Summary  of  boundary  conditions 


The  wall  conditions  are: 


6fj  =  0 


(92a) 


5ux  =  0 


(92b) 


63l  -  632  =  0 


(92c) 


(a)  for  the  transformed  wall  shear  method 
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£ 

as 


vv£ 


s-T-. 


6G  =  0 


(92d) 


(b)  for  the  skin  friction  coefficient  method 


(1  -  c1)b1  6G1  -  A1631  =  0, 


(92e) 


(c)  for  the  hybrid  method 

* 


(1  -  ex)  6G:  +  X2  =  P 

(d)  and  for  the  displacement  thickness  or  direct  method 

.2 


( 92  f ) 


6f2  -  &fl  -  p26u2  +  p1<5u1  +  ^ 


,2  (1  -  e2)  6G2  -  (1  -  ex)  «GX 


fl  "  f2  P1U1  +  P2U2  ’  12 


*  * 

b2  bl 

<2  ^2  '2  ^1 

n2  *  n  *  1 


"1 

1 


The  far  field  boundary  conditions  are: 

6uk+i  ■  0 

and,  from  the  two-point  spline  conditions, 

B116fK  +  B126UJ,  +  B13<5Gk  +  B14SPK 

+  AU6fK+1  +  A126uK+1  +  A135Gr+1  +  A14(5PK+1  =  T^, 


B316fK  +  B326UJ,  +  B33<5GK  +  B346^K 


+  A316fK+1  +  A325uk+1  +  A335GK+1  +  A34$PK+1  -  PK+1, 


B4L5fj,  +  B42<5uj,  +  B435Gj,  +  B446f3j, 

*  MI5fK+l  +  A«26u,.+1  ♦  A43«Gk+1  +  AMJf!^  -  0^ . 
For  the  displacement  thickness  method,  Eq.  (93b)  is  replaced  by 


6f 


K+l  +  A160K+1  =  nK+l  ‘  fK+l  ’ 

V2E, 


(5  /Re). 


w 


(92g) 


(93a) 


(93b) 


(93c) 


(93d) 


(93e) 
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For  the  direct  method,  the  final  boundary  condition  reduces  to 
5f3K+l  =  °- 


(93f) 


2.10  Block  Tridiagonal  Solver 

The  linearized  interior  spline  equations,  Eqs.  (79)  and  (81), 
together  with  the  wall  and  far  field  boundary  conditions,  Eqs.  (92)  and 
(93),  can  be  put  into  the  following  block  matrix  form: 


1Z1  +  C1Z2 

=  R  1 

i.Z.  .  +  A.Z.  +  C.Z.  . 

J  J’i  J  J  J  J+1 

=  Rj,  2  <  j  <  K, 

(94) 

BK+1  ZK  +  ^Sc+1  ZK+1 

=  ^4-1  ’ 

with  the  solution  vector 


<5f 

6u 

Z  *  6G 

j  53 

_  _  j 

and  the  right-hand-side  vector 


T 

L 

R  -  P 
j  Q 
_ j 

This  is  a  A  X  4  block  tridiagonal  system  which  can  be  accurately  and 
efficiently  solved  by  lower-upper  (L-U)  decomposition.  A  block 
tridiagonal  solver,  using  subroutines  developed  by  Blottner  [34]  which 
perform  partial  pivoting,  is  used  to  solve  the  matrix  equation  for  the 
correction  vectors  at  each  iteration.  Partial  pivoting  is  employed  to 
prevent  the  buildup  of  roundoff  errors. 


The  tridiagonal  system  of  Eq.  (94)  can  be  written  as 

AS  ■  5 

where  /a  can  be  decomposed  as 

A-1U 

with 

A  -  tBj,  A..  C.) 

IL  ■  iSj.  «j,  01 

1‘  -  to.  «j.  Vjl- 

The  forward  and  backward  substitution  equations  are  then 

Li  -  * 
rj  z  -  2 

Since  the  factorization  is  not  uniquely  determined  by  the  decomposition 
as  given,  a.j  in  ^  is  set  to  the  identity  matrix  (This  is  Keller  s  case 
(ii).  See  Ref.  [35]).  The  resulting  recursion  formulas  are 
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(95a) 

-v 

2  <  j  <  K+l 

W  ei  Vi- 

(95b) 

(95c) 

<5j  =  Cj  ,  1  <  j  <  K 

(95d) 

forward  and  backward  substitution  formulas  are 

Vi  ’  Ri' 

(96a) 

<5  .y .  -  R.  -  (3 .  y.  , ,  2  <  j  <  K+l 

J  J  J  J  J-1 

(96b) 

^+1  “  yK+l  * 

(96c) 

zj  ‘  yj  '  yj  zj-i-  "  -  J  -  1 

(96d) 

The  block  tridiagonal  solution  process  is  repeated  at  a  streamwise 
station  i,  updating  the  solution  variables  with  each  iteration,  until 


wamii  r.\ 


convergence  for  the  vector  of  solution  corrections  is  achieved.  The 
convergence  criterion,  applied  to  each  correction  at  each  q  nodal  point 
for  station  i,  is 

|g(n+i;)  .  gCn)  j  „  j5g(n)|  <  10'8f  l  <  j  <  K+l  . 

Upon  convergence  of  the  solution  at  a  streamvise  station,  the  iteration 


process  is  repeated  at  the  next  streamwise  station,  and  the  solution  is 
marched  progressively  downstream. 


2.11  Starting  Solution 

The  transformed  pseudo-self -similar  form  of  the  equations  allows 
the  solution  procedure  to  be  self-starting  since,  in  most  cases,  the 
solution  can  be  started  from  either  a  flat  plate  or  stagnation  point 
flow.  Both  of  these  flows  are  self-similar  in  nature,  and  the 
transformed  governing  equations  reduce  to  ordinary  differential 
equations  with  no  dependence  on  £. 

As  an  initial  guess  for  the  starting  solution,  a  fourth-order 
Pohlhausen  polynomial  is  used  to  approximate  the  self-similar  profiles 
of  (f,u,G,3)ij.  The  Pohlhausen  polynomial  has  the  form 

—  =  f  =  b£  +'cC2  +  dC3  +  eC4  (97) 

u  n 

e 

where  £  *  q/r^,  The  constants  b,c,d,  and  e  are  found  by  applying  the 

boundary  conditions  and  the  Falkner-Skan  equation 

f  +  ff  +  B ( 1  -  i2)  =  0  .  (98) 

qqq  qq  q 

The  polynomials  for  the  other  required  profiles  can  be  found  by 
integrating  Eq.  (97),  for  f,  and  differentiating  Eq.  (97),  for  G. 

An  initial  ^uess  for  3  can  be  found  by  substituting  the  previous 
results  into  the  Falkner-Skan  differential  equation,  yielding 


.’A  Li 
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3  -  ^  (b  -  2) 


(99a) 


1 1 

b  =  n  f 

<*>  w 


(99b) 


This  result  can  be  reversed  to  find  f  ,  with  (5  given,  for  the  direct 

w 

method  so  that 


t  i 

f  =  b/n_ 


( 100a) 


n_  P  +  12 


(100b) 


With  the  initial  conditions  provided  by  the  Pohlhausen  polynomial, 
the  starting  solution  is  obtained  by  solving  the  ODE  system  with  the 
standard  block  tridiagonal  method.  This  procedure  can  also  be  used  to 
obtain  any  Falkner-Skan  self-similar  solution. 

If  a  non-similar  solution  is  to  be  obtained,  the  starting  solution 
is  found  using  the  fourth-order  polynomial  approximation,  followed  by 
marching  the  solution  in  the  downstream  direction.  For  the  second 
streamwise  station,  two-point  backward  differences,  instead  of  the  usual 
three-point  differences,  are  used  to  approximate  the  £ -derivatives. 

Once  past  the  second  streamwise  station,  the  three-point  differences  are 
ordinarily  used. 


2.12  Determination  of  the  Edge  Velocity 

The  boundary- layer  edge  velocity,  ue,  is  required  to  obtain  the 
eddy  viscosity  distribution.  It  is  also  an  important  parameter  of  the 
boundary- layer  flow  and,  should  the  present  method  be  extended  to  a 
viscous-inviscid  procedure,  would  be  necessary  for  the  outer  iteration. 
For  these  reasons,  an  accurate  method  of  determining  the  edge  velocity 
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is  most  important.  The  edge  velocity  cannot,  however,  be  determined 
directly  from  the  solution  procedure,  since  the  pseudo-self-similar 
transformation  eliminates  ue  from  appearing  explicitly  in  the  governing 
equations.  For  this  formulation,  ue  must  be  obtained  through  an 
integration  of  the  solution  variable. 

The  pressure  gradient  parameter  is  defined  as 


__  du 
a  =  l£  __e 

P  u  d£ 
e 


(13) 


Noting  that  this  is  a  double  logarithmic  derivative,  Eq.  <  13)  can  be 
recast  as 


.  d(inu  ) 

1  o  =  _ e_ 

2  p  d(  in  C) 


(101) 


Rearranging  and  integrating  over  the  numerical  interval  i-1  to  i,  Eq. 
(101)  becomes 


in 


e. 

x 


1 


u 


ei-l 


r  /  0din£ 
l-l 


The  right-hand-side  of  the  previous  equation  may  be  approximated  using 
the  trapezoidal  integration  rule.  The  edge  velocity  at  station  i  is 
then  obtained  from 

r 


e . 
i 


e . 


(102a) 


i-1 


where 


.  (  n  ) 


r  *  |  (S(n)  +  3.  ,)  [nr 

A  i  i  - 1 


f  i  o ;  b ) 


i-1 


As  the  iteration  process  proceeds,  Eq  .  <  I'jZ)  must  be  updated  w ;  *  h  new 

values  of  3^  and 


5 


2 


* 


V*. 


i 


a 


a 


:  .1 

3 


vy 


1 


A 


<1 
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This  method  of  determining  ug  presents  a  problem  when  solving  a 
npnsimilar  flow  beginning  with  a  stagnation  point.  At  the  stagnation 
point,  i=l,  the  edge  velocity  is  zero.  To  calculate  subsequent  edge 


velocity  values,  the  value  of  ue  at  i=2  must  be  input  directly  or 
estimated.  This  anamolous  behavior  at  the  second  streamwise  station  is 
due  to  the  double  logarithmic  form  of  the  pressure  gradient  parameter 
definition. 


The  reason  for  this  behavior  is  apparent  if  the  flow  near  the 


stagnation  point  is  assumed  to  behave  like  a  wedge  flow  [36].  Then,  the 


edge  velocity  is  an  exponential  function  of  x  such  that  (3  =*  1) 


2-3 

u  =  cx  =*  cx 
e 


The  constant  c  can  not  be  determined  without  specifying  the  body  shape 


because  of  the  psuedo-self-similar  character  of  the  flow.  A  similar 
observation  is  made  by  Bradshaw,  Cebeci,  and  Whitelaw,  [37],  who  note 
that  it  is  necessary  to  specify  the  slope  c  near  a  stagnation  point. 


2.17  Normal  Coordinate  Distributions 

Although  the  governing  equations  are  formulated  to  allow  the  use  of 
the  adaptive  normal  coordinate,  the  use  of  this  coordinate  is  not  always 
practical  for  all  flows,  particularly  flows  which  are  entirely  laminar. 


For  laminar  flows,  the  system  is  returned  to  (£,  n)  coordinates  and  a 
more  straightforward  distribution  of  the  normal  coordinate  is  used. 

For  entirely  laminar  flows,  either  a  constant  n  stepsize  or  a 
geometric  progression  in  An  is  computed.  The  geometric  progression  is 


ib'-imed  from 
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Ar^  =  Ce  Anj_i»  2  <  j  <  K+l 
with  the  constant  ce  set  to  1.1. 

As  an  alternative  to  the  adaptive  normal  coordinate  for  turbulent 
boundary- layer  computations,  hyperbolic  stretching  functions  are  also 
installed  in  the  computer  code.  This  type  of  coordinate  distribution  is 
discussed  at  length  by  Vinokur  [38]  and  Thompson,  Warsi,  and  Mastin 
[39].  Here,  one-sided  hyperbolic  tangent  or  hyperbolic  sine  functions 
are  used  to  determine  a  coordinate  distribution  in  n.  The  details  of 
these  stretching  functions  are  not  included  here  since  they  are 
adequately  covered  in  the  noted  references . 


.  s*  ■ 


ti 
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Chapter  3 
RESULTS 


3 . 1  Laminar  Results 

The  following  subsections  describe  the  seven  laminar  test  cases 
which  were  run,  including  both  inverse  and  direct  solutions.  The  cases 
are  summarized  in  Table  2,  which  displays  pertinent  information  about 
each  run.  Geometric  progressions  were  used  for  the  normal  coordinate 
distributions  unless  otherwise  noted. 

3.1.1  Linearly  Decreasing  Mall  Shear 

The  linearly  decreasing  wall  shear  cases  were  originally  devised  by 
Keller  and  Cebeci  [8],  and  were  later  repeated  by  Cebeci  and  Keller 
[22],  Horton  [10],  and  Kaufman  and  Hoffman  [23].  The  cases  listed  here 
are  a  review  of  those  given  in  Reference  23  with  several  additions. 

The  cases  were  devised  to  specify  a  linearly  decreasing  transformed 
wall  shear  with  one  case  beginning  with  a  flat  plate, 

SiO  -  0.4696  (1  -  5).  5  >  0  , 
and  the  second  case  beginning  with  a  stagnation  point, 

S(0  =  1.23259  (1  -  C),  C  >  0  , 

with  both  proceeding  to  separation  at  E,  =  1.  For  this  case,  comparisons 
are  also  made  with  the  hybrid  solution  method,  the  displacement 
thickness  method,  and  a  direct  solution. 

Figure  1  gives  a  comparison  of  run  (la)  and  the  computational 
results  of  Keller  and  Cebeci  [8]  and  Cebeci  and  Keller  [22].  This 
figure  clearly  shows  the  breakdown  of  the  Cebeci  and  Keller  mechul 
function  near  separation.  The  results  from  the  present  method  compare 
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Table  2 

Boundary  Layer  Test  Cases 


Run  No. 

Mode 

B.D. 

M 

K 

Hoc 

Specifications 

Case  1:  Linearly  Increasing  Wall  Shear 

la 

I-A 

3  pt . 

23 

20 

7 

Flat  Plate 

lb 

I-A 

3  pt. 

23 

20 

7 

Stag.  Point 

lc 

I-A 

3  pt. 

23 

10 

7 

F.P. 

Id 

I-A 

3  pt. 

23 

20 

7 

F.P. ,  Const. An 

le 

D 

3  pt. 

23 

20 

7 

F.P. 

If 

I-G 

3  pt. 

23 

20 

7 

F.P. 

lg 

I-G 

3  pt. 

23 

20 

7 

S.P. 

lh 

I-E 

3  pt. 

23 

20 

7 

F.P. 

li 

I-E 

3  pt. 

23 

20 

7 

■ 

S.P. 

Case  2: 

Howarth  Flow 

2a 

D 

3  pt. 

22 

20 

7 

.... 

2b 

I-B 

3  pt . 

20 

20 

7 

— 

2c 

I-F 

3  pt . 

20 

30 

8 

.... 

Case  3: 

Potential  and  Viscous  Circular  Cylinder 

3a 

D 

3  pt. 

50 

I 

20 

i 

Potential 

3b 

I-A 

3  pt. 

4  7 

20 

7 

Potential 

3c 

D 

3  pt. 

40 

20 

7 

Hiemenz  Fit 

3d 

I-A 

3  pt. 

36 

20 

7 

Hiemenz  Fit 
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VI 

>*r 


VL 

>>| 

•a 

v| 


'.•I 


Run  No. 

Mode 

B.D. 

M 

K 

■ 


Specifications 


Case  4:  Horton  Parabolic  Wall  Shear  Distribution 


4a 

!  I-A 

3  pt. 

50 

30 

4b 

I-B 

3  pt. 

50 

30 

4c 

I-C 

3  pt. 

50 

30 

4d 

I-F 

2  pt. 

50 

30 

4e 

I-E 

2  pt. 

50 

30 

Case  5:  Analytic  6" /Re  Distribution 


5a 

5b 

5c 

5d 

5e 

5f 

5g 


Case  6 : 


I-F 

3  pt. 

50 

30 

I-F 

3  pt. 

50 

30 

I-F 

3/2  pt. 

50 

30 

I-F 

2  pt. 

50 

30 

I-C 

3/2  pt. 

50 

30 

I-F 

3  pt. 

100 

30 

I-F 

3/2  pt. 

100 

30 

_ _ 

Horton/Ntim  5  /Re 

Distribution 

D 

3  pt. 

76 

30 

I-F 

3  pt. 

7  6 

30 

r-F 

3/2  pt. 

76 

30 

6  *  5.6 

max 
* 

6  =8.6 

max 
A 

6  =  8.6 

max 

* 

6  =  8.6 

max 

* 

6  =  5.6 

max 

* 

5  =  5.6 

max 

a 

6  =8.6 

max 
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Table  2  (Continued) 


Run  No. 

Mode 

B.D. 

M 

K 

H 

Specif ications 

Case  7:  Klineberg-Steger  t 

Distribution 

7a 

I-D 

3  pt. 

120 

20 

7 

sep.  2-6 

7b 

I-D 

3  pt. 

120 

20 

7 

sep.  2. 5-5.5 

7c 

I-D 

3  pt. 

120 

20 

7 

sep.  3-5 

Modes :  D-Direct  , , 

I-Inverse  A  -  f  specified 

B  -  C^  /Re  specified 
e 

C  -  C^  /Re  specified 

oo 

D  -  t  specified 

E  -  6  specified 
* 

F  -  6  /Re  specified 


G  -  Hybrid  Method 
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well  with  those  of  the  Keller  and  Cebeci  eigenvalue  method  to  5  =  0.9, 
and  proceed  through  separation  while  the  other  two  methods  fail. 

Figure  2  gives  a  comparison  of  results  for  the  transformed  wall 
shear  method,  the  hybrid  method,  and  the  displacement  thickness  method. 
All  agree  quite  well  up  to  separation,  where  the  hybrid  method  fails  to 
converge.  The  other  methods  continue  through  separation  with  excellent 
agreement.  It  should  be  noted  that  calculations  are  only  shown  slightly 
past  the  separation  point  since  the  flow  separation  is  severe  and  the 
computation  methods  fail  for  large  regions  of  backflow. 

A  study  of  the  accuracy  of  the  fourth-order  splines  for  different 
normal  coordinate  distributions  is  displayed  in  Figure  3.  Here,  a 
comparison  is  made  between  a  twenty-point  geometric  progression  in  n,  a 
twenty-point  constant  increment  in  n»  and  a  ten  point  geometric 
progression  in  n .  All  three  resulting  8  distributions  are  identical, 
indicating  the  ability  of  the  splines  to  yield  accurate  results  with 
significantly  fewer  nodal  points. 

The  direct  method  was  tested  on  the  fiat  plate  to  separation 
problem  by  using  the  ue  distribution  obtained  from  run  (la).  The 
transformed  wall  shear  from  this  run  is  compared  in  Figure  4  with  direct 
results  for  Reference  8  as  well  as  with  the  exact  value  specified  in  run 
(la).  The  comparison  is  quite  good  with  only  slight  deviation  of  the 
present  results  from  the  exact  distribution  at  separation,  which  is 
expected  for  a  direct  boundary- layer  method. 

The  calculations  of  the  various  inverse  methods  for  the  stagnation 
point  to  separation  case  are  shown  in  Figure  5,  with  comparison  of 
computational  results  from  Reference  8.  Again,  the  agreement  is  quite 


i-igure  4.  Comparison  of  fw  from  Direc 
Case  1.  flat  Plate  to  Separ 


*4 
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good,  with  both  the  wall  shear  and  displacement  thickness  methods 
continuing  through  separation. 

3.1.2  Howarth  Flow 

A  boundary  layer  with  a  linearly  decreasing  edge  velocity,  the 
Howarth  flow  [1,40],  has  become  a  well-known  test  case  for  boundarv- 
layer  computations,  and  is  included  here  mainly  as  a  test  for  the  direct 
method  since  no  separation  data  is  available.  Two  inverse  runs  are 
included  for  completeness. 

Figure  6  shows  both  direct  and  inverse  methods  to  have  excellent 
agreement  with  both  new  (Edwards  et  al.  [29])  and  old  (Smith  and  Clutter 
[1])  computational  methods.  A  tabular  comparison  of  dimensionless  wall 
shear  values  obtained  by  the  present  method  and  results  from  Howarth 
[41],  Cebeci  and  Smith  [33],  and  Smith  and  Clutter  [1]  is  presented  in 
Table  3.  Figures  7  and  8  compare  the  exact  edge  velocity  and  pressure 
gradient  parameter  distributions,  respectively,  with  the  inverse 
results.  The  agreement  is  excellent,  with  only  slight  deviation  evident 
in  S  near  separation.  This  discrepancy  is  most  likely  due  to  the  direct 
solution  being  used  as  input  for  the  inverse  runs. 

3.1.3  Circular  Cylinder 

As  for  the  Howarth  flow,  the  flow  past  a  circular  cylinder  has  also 
become  a  standard  test  case  for  boundary  layer  methods.  Here  direct 
solutions  are  computed  using  the  edge  velocity  distribution  from 
potential  theory  [40]  as  well  as  using  the  Hiemenz  polynomial  fit  of  the 
edge  velocity  for  viscous  flow.  This  case  is  included  to  test  the 
ability  of  the  direct  and  inverse  methods  to  obtain  a  solution  beginning 
with  a  stagn>tion  point.  Although  Case  1  included  a  flow  beginning 
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with  a  stagnation  point,  the  linearly  decreasing  wail  snear  problem  is 
artificial,  created  only  to  test  a  solution  procedure,  and  does  not  have 
a  unique  edge  velocity  distribution. 

Figure  9  gives  a  comparison  of  transformed  wa i 1  shear  obtained  hv 
the  present  direct  solution  method  with  Terrill  (42]  for  the  potential 
flow,  and  Smith  and  Clutter  (Li  for  the  Hiemen.:  poivnomal  f  i  r  .  The 
agreement  is  quite  acceptable  for  both  cases,  with  the  present  mefh  : 

I  » 

yielding  slightly  larger  values  of  f  than  References  i  and  *2.  A 
comparison  of  v'Re  distributions  for  the  potential  flow  case  is  give-, 

<30 

in  Figure  10.  The  values  obtained  by  the  direct  and  inverse  run,  <  a 
and  (3b)  compare  very  wed  with  results  given  by  Cebeci  and  3r -ids haw 
[43],  The  direct  solutions  predict  separation  at  approx  mate iv  10  3’, 
for  the  potential  f  ,  and  80°,  for  the  viscous  flow.  These  values  are 
generally  accepted  as  the  correct  locations  for  separation  [40';. 

The  results  of  the  inverse  calculations  are  compared  with  the  exact 
ue  distribution  as  given  by  White  [40],  and  the  anaivt ieai Iv 

t  * 

determined  5  distribution  in  Figures  11  and  12.  Here,  f  distributions 
obtained  from  the  direct  solutions  ( 3a )  and  (3c)  are  used  as  input  f >r 
the  inverse  solution  procedure.  The  match  between  the  inverse  solutions 
and  the  analytical  curves  is  almost  exact. 

3.1.4  Horton's  Parabolic  Wall  Shear  Distribution 

f  « 

Horten  ( 1 0  J  devised  a  parabolic  distribution  of  f^  to  test  the 
ability  of  the  boundary-layer  equations  to  treat  small  regions  of 
separation.  The  distribution  of  transformed  wall  shear  is  given  by 
SCO  =  0.4696  (1  -  O  (1  -  0.526490,  £  >  0. 

Separation  occurs  at  (  *  1,  with  reattachment  at  £  *  1.9.  Because  of 
its  simplicity  and  the  existence  of  a  small  separation  bubble,  this  case 
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was  used  as  a  test  for  all  five  o_  the  inverse  options.  The  transformed 
wall  shear  run,  (4a),  was  used  to  generate  input  data  for  the  other  four 
runs . 

Figure  L3  compares  Horton's  computational  data  with  three  inverse 
"  , —  *  . — 

options,  specifying  fy  ,  Cf  /Re,  and  <5  v^Re,  from  the  present  study. 

e 

The  agreement  among  all  results  is  very  good,  although  the  displacement 
thickness  run  displays  a  slightly  different  trend  from  the  other  cases. 

Two  interesting  discoveries  were  made  while  obtaining  the 
displacement  thickness  runs.  Accurate  results  could  not  be  obtained 
without  overestimating  the  value  of  rjoo  used  for  the  wail  shear  cases  by 
ten  to  fifteen  percent,  and  using  two-point  backward  differences, 
instead  of  the  usual  three-point  differences,  over  the  ^-interval  where 
the  flow  separates.  Oscillations  in  the  normal  velocity,  V,  occurred 
when  the  case  was  run  using  three-point  backward  differences  for  the 
entire  domain.  The  cases  where  wall  shear  or  skin  friction  coefficient 
distributions  were  specified  did  not  exhibit  these  oscillations. 

Figure  14  compares  runs  (4b)  and  (4c),  where  /Re  and  /Re 

0  oo 

are  specified,  respectively.  Figure  15  makes  a  similar  comparison 
between  the  displacement  thickness  cases.  For  both  figures,  the 
compared  curves  are  identical,  with  the  exception  of  a  small  oscillation 
for  the  Cj.  /Re  run  near  the  starting  solution.  This  oscillation  does 

OO 

not  affect  the  solution  downstream,  and  its  origin  is  unclear. 

Figure  16  compares  the  computed  5  /Re  distribution  from  run 

(4a)  with  the  computational  data  of  Horton  [10].  The  agreement  is 

exceptionally  good,  with  essentially  no  deviation  between  the  two 


curves . 
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Figure  17.  Comparison  of  fw  from  Inverse  Methods  with  Exact 
Distribution  for  Case  4. 
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I  I 

Computed  values  of  f  are  compared  in  Figures  17,  18,  and  19.  In 
Figure  17,  the  skin  friction  coefficient  and  displacement  thickness 
methods  are  compared.  Again  the  curves  show  minor  but  acceptable 
differences  between  the  methods,  particularly  in  the  latter  half  of  the 
domain.  The  results  of  the  skin  friction  coefficient  method,  however, 

I  t 

exactly  match  the  exact  values  of  f  as  specified  in  run  (4a).  Figure 
18  compares  the  skin  friction  coefficient  runs,  with  no  deviation 
between  results,  ignoring  the  small  fluctuation  near  the  starting 
solution.  Figure  19  shows  similar  agreement  for  the  displacement 
thickness  runs. 


3.1.5  Analytic  Displacement  Thickness  Distribution 

Carter  [13]  developed  an  algebraic  distribution  of  the  untransformed 

*  / — 

displacement  thickness  6  vRe,  which  matches  the  flat  plate  displacement 
thickness  at  the  start  and  end  of  the  distribution  and  reaches  a  maximum 
value  at  a  specified  interior  location.  Carter  uses  this  distribution 
to  test  a  forward  matching  scheme  with  FLARE  as  well  as  a  global 
iteration  scheme.  Cebeci,  Keller,  and  Williams  [7]  use  the  same 
distribution  to  test  the  marching  solution  with  the  iterative  DUIT 
procedure  for  separation  bubbles.  The  distribution  as  given  by  Cebeci 
et  al.  [7]  is  (with  errors  corrected): 


1 . 7208(x) 


1  £  x  £  xi 


6  (x)  =  »/Re  6  (x)  =  (a,  +  a~(x-x,)  +  an(x-x,)Z  +  a,(x-x.)  ,  x,  <  x  < 

I  1  Z  1  31  41  1—  — 


a1  +  a2(x-x2)2  +  a2(x-x2)3, 


x2  £  x  £  x3 


where 


a3  =  1.7208(x1) 


a2  =  (0.5)  (1.7208)/(x1) 


88 


3u  *  — r  [“  (6  -  a  )  -  4a  ] 
3  (24^)  4^  max  o  1 


2  ,a2  A1  ,t*  x  , 

a.  *  — r~  [ — r - (6  -  a.  )  } 

4  (a3}  2  max  1 


-A 

a,  =  6 
1  max 


a7  =  -  <-f)  [3(<S*  -  2.25)] 

4  .  Z  max 

*2 


“3  *  (-f  >  '2<Cx  -  2'25)1 
A2 


—  X2  _  ^2  *  x3  ”  x2 


x^  =  1.065,  =  1.35,  =  1.884 


Figure  20  displays  the  resulting  distributions  for  ^max  *5.6  and 

—it 

6  *  8.6. 

max 


Figure  21  gives  a  comparison  between  the  input  displacement 
thickness  distribution  and  the  generated  distribution  from  run  (5a), 
specifying  6  ,  and  run  (5e),  specifying  >/Re.  The  agreement  between 

QO 

the  displacement  thickness  solution  and  the  exact  distribution  is  quite 
good  overall  but  overshoots  where  the  distribution  matches  the  flat  plate 
value  downstream.  The  skin  friction  case  provides  correct  results  while 
the  flow  remains  attached,  but  degenerates  rapidly  upon  separation  and 
fails  inside  the  separation  bubble.  The  cause  of  this  failure,  in  light 
of  the  apparent  successful  calculation  for  the  Horton  parabolic 
distribution  [10],  will  be  discussed  in  the  following  chapter. 

The  skin  friction  distribution  resulting  from  run  (5a)  is  compared 
with  computational  results  of  Carter  [13]  and  Cebeci  et  al.  [7]  in  Figure 
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22.  The  agreement  within  the  region  of  separation  is  particularly  good, 
with  the  present  method  matching  the  global  iteration  results  of  Carter 
and  the  DUIT  procedure  of  Cebeci  et  al.  more  closely  than  the  forward 
marching  results  of  Carter.  Figure  23  compares  the  edge  velocity 
distribution  for  the  same  case  with  the  global  iteration  results  of 
Carter.  Again,  the  comparison  is  quite  good,  although  the  region  of 
overshoot  is  clearly  identifiable  following  reattachment. 

Figure  24  displays  two  computational  runs,  (5b)  and  (5c)  for  the 

analytic  6  distribution  where  <5  =8.6.  Run  (5b)  was  made  with 

max 

three-point  backward  differences,  as  was  run  (5a).  For  the  larger 
displacement  thickness  case,  however,  it  was  necessary  to  switch  to  two- 
point  backward  differences  upon  separation  to  obtain  an  accurate 
solution.  Using  two-point  backward  differences,  the  solution  matches 
the  analytic  6  distribution  very  well,  although  the  overshoot  behavior 
is  still  apparent  past  x=1.7.  For  this  case,  the  overshoot  causes  the 
solution  to  diverge  near  x=1.8.  This  behavior  following  reattachment  is 
discussed  later  in  this  section  and  in  Chapter  4. 

A  comparison  of  C^  /Re  distributions  is  made  in  Figure  25.  From 

oo 

the  present  study,  we  compare  displacement  thickness  method  results  with 
three-point,  three/ two-point  combination,  and  two-point  backward 
differences.  The  three-point  backward  difference  result  diverges  rapidly 
following  separation.  The  three/two-point  combination  and  the  two-point 
backward  difference  runs  are  essentially  identical.  These  results  are 
compared  with  forward  marching  and  global  iteration  results  from  Carter 
[13]  and  forward  marching  results  with  DUIT  from  Cebeci  et  al.  [7). 

Again,  the  present  results  match  the  global  iteration  or  DUIT  results 
more  closely  than  the  forward  marching  results  with  FLARE  provided  by 


97 


Carter.  The  overshoot  behavior  is  also  more  evident  here  than  in  run 
(5a).  The  edge  velocity  distribution  from  the  same  runs  is  compared  with 
Carter's  global  iteration  results  in  Figure  26.  Here,  the  calculations 
of  the  present  method  fall  slightly  below  the  results  of  Carter,  but  only 
by  approximately  one  percent.  The  overshoot  behavior  is  very  apparent 
downstream  of  x=1.6,  shifting  the  minimum  value  of  ue  upstream 
approximately  0.05  units. 

The  suspected  cause  of  the  overshoot  behavior  is  streamwise 
discretization  error.  To  study  this  effect,  the  number  of  streamwise 
points  was  doubled,  and  the  cases  rerun.  Figures  27  and  28  display  the 
results  for  the  two  maximum  displacement  thickness  values.  In  both 
cases,  and  particularly  for  =  8.6,  there  is  a  marked  improvement 

over  the  original  results. 

Figure  29  shows  the  remarkable  improvement  obtained  with  the 
decrease  in  streamwise  stepsize  for  the  skin  friction  results.  The 
results  of  run  (5g)  match  the  global  iteration  and  DUIT  results  much  more 
closely  than  before,  with  only  minor  overshoot  evident  compared  to  the 
initial  run  with  50  streamwise  grid  points.  Figures  30  and  31  show 
similar  results  for  the  edge  velocity  distribution,  again  with  the  most 
improvement  shown  for  the  case  where  =  8.6.  Although  there  is 

still  a  slight  upstream  shift  near  reattachment  evident  in  Figure  31,  the 
decrease  in  stepsize  dramatically  effected  the  results  throughout  the 
domain. 

A  series  of  plots  are  also  included  here  detailing  the  streamwise 
velocity  profile  as  the  solution  marches  downstream  through  separation 
for  run  (5a).  Figures  32  a-c  show  velocity  profiles  from  the  separation 


Figure  28.  Effect  of  Streamwise  Stepsize  Reduction.  Comparison  ot 
6*  /Re  for  Case  5  (5*  *  8.6). 
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Figure  29.  Effect  of  Streamwise  Stepsize  Reduction.  Comparison 

of  Cc  >/Re  for  Case  5  (5*  *  8.6). 
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Separation  Sequence.  Streamwise  Velocity  Profiles  at 

Three  Streanwise  Stations  for  Case  5  (6*  =5  63 
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(a)  Station  7 


RUN  (5a) 
STATION  30 
X  =  1.54002 


Figure  33. 


Wall  Region  Detail  of  Separation  Sequence.  Streamwise 


Velocity  Distribution  at  Three  Streamwise  Stations  for 
Case  5  (6*ax  =  5  ■  6)  • 


(a)  Station  7 


Wall  Region  Detail  of  Separation  Sequence.  Streacvise 
Velocity  Distribution  at  Three  Streamwise  Stations  for 


Case  5  (c*ax  =  5.6) . 


(b)  Station  15 


Figure  33.  Wall  Region  Detail  of  Separation  Sequence.  Streamwise 

Velocity  Distribution  at  Three  Streamwise  Stations  for 

Case  5  (6*  =  5.6). 
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(c)  Station  30 
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point  to  the  maximum  negative  value  of  f  .  Figures  33  a-c  provide 
detailed  views  of  the  backflow  regions  of  the  same  profiles. 

3.1.6  Horton  Displacement  Thickness  Distribution 

An  inverse  case  using  the  experimental  data  of  Ntim,  as  given  by 
Horton  [11],  and  the  computational  data  of  Horton  from  the  same 
reference,  was  attempted  using  the  displacement  thickness  method.  The 
original  data  gathered  by  Ntim  was  obtained  by  placing  a  circular 
cylinder  over  a  flat  plate  to  alter  the  pressure  distribution  and  cause 
separation.  Following  Horton,  the  direct  method  was  used  to  generate  a 
displacement  thickness  distribution  to  a  location  near  separation.  The 
experimental  <5  /Re  distribution  of  Ntim  was  then  fitted  to  the 
computational  distribution  so  that  a  complete  inverse  solution 

vU  _ 

could  be  obtained.  Figure  34  shows  the  6  /Re  distribution,  resulting 
from  the  direct  solution  and  the  matched  experimental/computational 
distribution,  used  as  input  to  the  inverse  method.  Figure  35  is  a  detail 
of  the  direct  distribution  and  the  composite  distribution  around  the 
matching  point.  The  combined  distributions  are  joined  using  a  cubic 
spline  interpolation  routine. 

Jw  _ 

Figure  36  displays  the  input  6  /Re  distribution  with  the  generated 
<5~  /Re  d  istributions  from  two  inverse  runs,  (6b)  and  (6c).  As  was 
discovered  in  previous  cases,  two-point  backward  differences  had  to  be 
used  once  the  flow  separated  to  obtain  a  solution  into  the  separated 
zone.  A  noteworthy  feature  is  the  deviations  between  the  input 
displacement  thickness  and  the  generated  distribution  near  the  matching 
point.  A  similar  result  is  apparent  in  the  edge  velocity  distribution, 
shown  in  Figure  37,  which  peaks  above  the  edge  velocity  distribution  used 
as  input  to  the  direct  solution. 
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The  cause  of  this  discrepancy  becomes  more  apparent  after 


considering  the  plot  of  skin  friction  coefficients  in  Figure  38.  In 


this  figure,  several  discontinuities  are  easily  visible  near  the  6  ^Re 


distribution  matching  points.  These  peaks  in  C,  v^Re  arise  because  only 


the  displacement  thickness  was  forced  to  be  continuous  at  the  matching 
points,  and  not  its  derivatives.  This  procedure  was  also  used  in 
obtaining  the  edge  velocity  distribution  for  the  direct  solution.  Some 
evidence  of  the  ug  distribution  matching  is  evident  in  the  skin  friction 
coefficients  generated  by  the  direct  solution.  The  inverse  solution, 
however,  appears  to  be  much  more  sensitive  to  such  matchings.  The 
mismatched  slopes  affect  the  solution  because  the  displacement  thickness 
is  related  to  the  streamf unction.  The  discontinuities  become  apparent 


in  the  skin  friction  coefficient  since  C,  is  a  function  of  f  ,  the 

f  w 


second  derivative  of  the  streamfunction,  and  differentiation  tends  to 
amplify  any  noise  in  the  calculation  procedure. 


3.1.7  Klineberg-Steger  Wall  Shear  Distribution 

An  analytic  distribution  of  wall  shear  was  devised  by  Klineberg  and 
Steger  [9]  to  test  an  inverse  scheme  which  they  developed.  The 


distribution  specifies  x,  which  is  related  to  fy  by 


_  /2£_ 


fw (5>  */r;  T(x)  • 
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The  algebraic  distribution  is  given  by 


'  0  332 

x(x)  =  fT -  (X“X1>(X-X2)  *  T  ^  ( X )  ,  0  <  X  <  Xp  X  >  Xp 


KS 


t(x)  =  x^xjfl  +  a(x-x1)(x-x2>  ] ,  <  x  <  x^. 
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Previous  calculations  show  that  the  wall  shear  methods  work  well  for  some 


flow  situations  (Horton's  parabolic  shear  distribution,  subsection  3.1.4) 
but  fail  completely  for  other  situations  where  they  are  expected  to  work 
(Carter's  analytic  displacement  thickness  cases,  subsection  3.1.5).  The 
Klineberg-Steger  distribution  was  used  in  this  investigation  to  provide  a 
better  qualitative  understanding  of  what  conditions  cause  the  wall  shear 
cases  to  fail. 

One  of  the  original  cases  run  by  Klineberg  and  Steger  with  the 
above  distribution  of  t  sets 
CRS  =  12.0,  a  =  0.0 
=  2.0,  x^  =  6.0 

For  this  investigation,  this  case  appears  as  run  (7a).  The  constants 

and  x^ ,  x^  were  then  changed  for  two  other  runs  so  that  the  slope 
d  '  ' 

—  (f^  )  into  the  region  of  separation  was  decreased.  This  also 
decreased  the  size  of  the  region  of  separation.  For  the  second  and  third 
computations,  the  constants  were  set  at: 

(1)  CKS  =  13.75,  xx  =  2.5,  x2  =  5.5 

(2)  CKS  -  15.0,  x,  -  3.0,  x2  =  5.0  . 

This  approach  was  taken  since  it  was  suspected  that  the  magnitude  of  the 

I  t 

f  slope  into  the  separation  region  was  a  main  cause  of  the  wall  shear 
inverse  option's  inability  to  successfully  treat  some  cases. 

The  three  cases  run  are  shown  in  Figures  39  and  40.  Figure  39 

f  » 

displays  the  distributions  of  f  .  Run  (7c),  which  has  the  smallest 

d  '  ' 

value  of  —  (f  )  at  separation,  passes  smoothly  through  separation  and 
reattaches  downstream.  This  case  has  the  smallest  separation  bubble, 
extending  from  x=3.0  to  x=5.0.  The  middle  case,  which  is  separated  from 
x=2.5  to  x=5.5,  fails  to  converge  in  the  separated  region  at  x=3.7, 


before  reaching  to  the  maximum  value  of  f  .  The  most  severe  case,  run 

w 

(7a),  fails  to  converge  properly  immediately  after  reaching  separation. 
The  pressure  gradient  solutions  for  the  three  cases  are  plotted  in  Figure 
40.  A  discussion  of  the  behavior  of  the  wall  shear  method  is  included  in 
the  following  chapter. 

3 . 2  Turbulent  Results 

The  application  of  the  present  method  to  turbulent  boundary  layers 
proved  to  be  disappointing.  As  a  test  case,  both  the  direct  and  inverse 
methods  were  us  d  to  compute  a  turbulent  boundary  layer  on  a  flat  plate. 

A  transition  intermittency  factor,  as  described  in  Appendix  B,  was  used 
to  permit  a  controlled  transition  to  a  fully  turbulent  flow  from  a 
laminar  starting  solution.  A  hyperbolic  tangent  stretching  function  was 
used  to  generate  a  normal  coordinate  distribution  for  the  flat  plate 
test . 

The  resulting  direct  solution  proved  to  be  very  unstable,  and 
diverged  only  several  steps  downstream  after  becoming  fully  turbulent. 

The  number  of  iterations  required  for  convergence  continued  to  increase 
as  the  solution  marched  downstream  even  after  the  full  value  of  the  eddy 
viscosity  was  reached  in  the  computations.  The  solution  was  permitted  up 
to  fifty  iterations  to  converge. 

The  attempted  inverse  solution  proved  to  be  more  unstable  than  the 
direct  solution,  and  diverged  only  a  short  distance  into  transition.  The 
skin  friction  distribution  generated  by  the  direct  method  was  used  as 
input  for  the  skin  friction  coefficient  solution  method.  Again,  a  large 
number  of  iterations  were  permitted. 

In  an  effort  to  locate  the  cause  of  the  divergent  behavior,  two 
other  similar  but  less  complicated  eddy  viscosity  models  were  coded  and 
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compared  with  the  Baldwin-Lomax  model.  The  standard  Cebeci-Smith  model 
[33],  in  basic  form  sans  corrections,  and  a  simple  model  applied  by 
Carter  and  Wornora  [14]  were  used.  Both  models  yielded  eddy-viscosity 
distributions  almost  identical  to  those  produced  by  the  Baldwin-Lomax 
model  for  the  flat  plate. 

An  attempt  was  also  made  at  increasing  the  transition  distance  by 
increasing  the  constant  in  the  transition  intermittency  factor.  This 
only  slightly  prolonged  the  convergence  of  the  solution  downstream. 

After  a  complete  check  of  all  analysis  and  program  coding  with  no 
errors  being  found,  the  attempt  at  the  turbulent  solution  was  abandoned. 
Future  work  will  be  needed  to  determine  the  cause  of  the  present  method's 
inability  to  handle  the  turbulent  case. 
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Chapter  4 


DISCUSSION 


A.l  Inverse  Methods 


In  the  present  investigation,  four  different  inverse  methods  for 


boundary- layer  flows  have  been  studied.  From  examining  the  results,  it 


is  clear  that  for  laminar  flows,  these  methods  work  very  well,  in 


general.  The  methods  in  their  current  formulations  are  not  without 


their  problems,  however,  and  for  turbulent  cases,  the  problems  are 


severe.  The  following  sections  discuss  the  numerical  aspects  of  the 


test  cases  given  in  the  previous  chapter  in  relation  to  the  various 


methods  and  formulations  used  in  this  investigation. 


The  inverse  methods  consist  of  the  transformed  wall  shear  method. 


the  skin  friction  coefficient  methods,  the  hybrid  method,  and  the 


displacement  thickness  method.  Each  of  these  methods  is  discussed 


individually  in  the  following  subsections. 


It  should  be  noted  that,  aside  from  Horton  [10,11],  this  study  is 


the  only  one  known  to  the  author  which  formulates  the  inverse  problem  in 


full  (streamwise  and  normal)  Levy-Lees  coordinates.  This  study  differs 


from  those  of  Horton  in  the  discretization  and  solution  procedure 


followed  once  the  equations  are  cast  in  pseudo-self-similar  form.  No 


other  studies  exist,  to  the  author's  knowledge,  which  cast  the  turbulent 


boundary- layer  equations  in  this  form  (modified  Levy-Lees)  for  an 


inverse  solution. 


£ 


12  2 

4.1.1  Wall  Shear  Methods 

The  transformed  wall  shear  method  proved  to  be  the  least 

complicated  of  the  inverse  formulations,  requiring  an  easy  to  apply 

additional  boundary  condition  at  the  wall.  For  all  cases  computed  using 

this  boundary  condition,  quadratic  convergence  was  observed,  with  only 

three  iterations  required  typically  to  satisfy  the  convergence  criteria. 

The  methods  allowing  the  specification  of  the  skin  friction 

coefficients  and  are  more  practical  from  a  physical  standpoint, 
e  00 

These  methods,  however,  require  the  development  of  the  complex  boundary 
conditions  described  in  Sections  2.4.2  and  2.9.3.  The  use  of  these 
boundary  conditions,  which  relate  the  quantities  £  and  ue  to  the 
pressure  gradient  parameter  correction,  is  necessary  to  accelerate  the 
convergence  of  the  solution.  With  these  boundary  conditions  installed, 
quadratic  convergence  is  maintained.  Without  the  correction  form  of 
these  boundary  conditions,  convergence  is  achieved  typically  only  after 
ten  or  more  iterations. 

The  wall  shear  methods  do  experience  a  problem  with  some  separated 
flows.  Flows  which  decelerate  rapidly  from  positive  to  significant 

i  i 

negative  values  of  f  over  a  small  streamwi se  distance  can  cause  the 
wall  shear  methods  to  fail.  Upon  examining  the  distributions  of 
streamfunction,  f,  and  velocity,  u,  across  the  boundary  layer,  an 
additional  reflex  in  the  stream  function  profile  is  found  near  the  wall, 
with  an  erroneous  positive  value  of  f  directly  adjacent  to  the  wall. 

This  eventually  leads  to  oscillations  in  the  solution.  A  comparison  of 
incorrect  profiles  oi  f  and  u  near  the  wall  and  correct  profiles  for  a 
similar  separated  region  are  shown  in  Figures  41  and  42.  This  failure 


Incorrect  Profiles  Near  Wall  from  Skin  Friction 
Coefficient  Option. 


(b)  Normalized  Velocity,  f' 
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is  apparently  inherent  in  the  spline  formulation  in  combination  with  the 
boundary  conditions  applied  at  the  wall  for  the  wall  shear  methods. 

Based  upon  a  comparison  of  results  from  the  wall  shear  and 
displacement  thickness  methods,  it  is  clear  that  there  is  a  difference 
in  formulation  which  must  account  for  the  displacement  thickness 
method's  ability  to  successfully  calculate  regions  of  separation  where 
the  wall  shear  method  fails.  This  difference  is  the  form  of  the  wall 
boundary  conditions  for  the  two  methods.  The  wall  boundary  conditions 
for  the  wall  shear  methods  only  involve  the  variables  f,  u,  G  at  j  =  1, 
with  no  dependence  on  j  =  2  (See  subsection  2.8.1).  Only  the  condition 
on  8  involves  values  at  both  j  =  1  and  j  =  2,  but  this  condition  has  no 


effect  on  the  profiles  for  f,  u  and  G. 

The  displacement  thickness  case  does  include  a  two-point  spline 
boundary  condition  at  the  wall  which  links  the  values  of  f,  u,  and  G  at 


points  i  =  1  and  i  =  2.  This  condition  apparently  forces  the  splines  to 
behave  in  the  correct  manner  near  the  wall  when  the  flow  separates, 
since  the  displacement  thickness  case  was  used  successfully  on  all 
occasions.  The  lack  of  this  two-point  condition  in  the  wall  shear 


methods  allows  the  splines  to  "wiggle"  near  the  wall.  This  phenomenon 

d  '  ' 

seems  to  occur  when  the  value  of  n— (f  )  exceeds  some  undetermined  limit 

dx  w 

near  or  following  separation.  This  conclusion  is  based  upon  the 
qualitative  study  of  laminar  case  number  7.  This  case  also  demonstrates 

I  I 

that  the  wall  shear  methods  fail  when  the  value  of  f  exceeds 

w 

approximately  -0.1  or  the  backflow  velocities  become  larger  than  one 
percent  of  the  edge  velocity.  The  inclusion  of  the  two-point  boundary 
condition  for  the  displacement  thickness  method  may  be  the  reason  this 

i  i 

method  can  also  handle  separated  zones  with  values  of  f  much  larger 


128 


than  the  wall  shear  methods  can  tolerate.  This  explains  why  the  wall 
shear  methods  perform  well  for  mild  cases,  such  as  the  Horton  parabolic 
distribution  (case  number  4)  but  fail  for  other  problems,  such  as  the 
separation  bubble  generated  by  Carter's  analytical  distribution  of 
displacement  thickness  (case  number  5). 

In  the  current  formulation,  however,  it  is  not  possible  to  install 
a  two-point  spline  boundary  condition  at  the  wall  for  the  wall  shear 
methods.  Considering  the  conditions  that  exist  at  the  wall,  only  the 
condition  on  8  may  be  transferred  to  the  outer  edge,  since  the 
homogeneous  conditions  on  f  and  u  must  remain  at  the  wall,  as  must  the 
inverse  wall  shear  condition  itself.  Unfortunately,  it  was  found  that 
transferring  the  two-point  condition  on  8  to  the  far  field  produced 
oscillations  there  which  occurred  for  attached  and  separated  flow 
conditions  and  which  significantly  decreased  marching  distance.  The 
conclusion  is  that,  in  the  present  formulation,  the  required  two-point 
wall  boundary  condition  which  would  allow  the  wall  shear  methods  to 
accurately  handle  larger  regions  of  separation  can  not  be  implemented. 

Despite  the  inability  of  the  wall  shear  methods  to  handle  larger 
regions  of  separation,  the  methods  have  proved  to  be  very  accurate  for 
milder  separated  regions  and  for  attached  flows.  This  point  is  well 
illustrated  by  the  results  for  cases  1-4,  where  the  wall  shear  methods 
provide  excellent  results,  particularly  for  Horton's  parabolic  shear 
distribution. 


4.1.2  Hvbrid  Method 

The  hybrid  method  described  in  Section  2.4.3  is  an  inverse  method, 
developed  solely  for  this  study,  which  takes  a  somewhat  different 
approach  to  specifying  5  than  other  displacement  thickness  methods, 
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such  as  those  presented  by  Carter  [12]  or  Cebeci  [15].  Initially  the 
inverse  procedure  presented  here  was  developed  as  a  wall  shear  method; 
and  hence  the  hybrid  boundary  condition  was  devised  to  make  use  of  the 
wall  shear  condition  for  specification  of  a  displacement  thickness. 

The  form  presented  in  this  study  has  not  been  evolved  to  the  same 
extent  as  the  wall  shear  and  displacement  thickness  methods,  since  it 
was  never  completely  linearized.  Thus  the  method  tends  to  require  more 
iterations  to  achieve  convergence  than  the  other  methods  discussed.  A 
standard  iteration  count  at  a  single  streamwise  station  is  typically 
between  five  and  twelve,  depending  on  the  condition  of  the  boundary 
layer  at  that  location. 

A  major  problem  with  the  hybrid  boundary  condition  is  the  saddle 

point  occurring  in  the  H-g  curve  near  separation.  The  value  of  the 
dH 

derivative  ^  for  mild  adverse  pressure  gradients  is  near  -1.0.  As  the 

dH 

saddle  point  is  approached  from  upstream,  however,  the  value  of  — 
approaches  negati'  a  infinity,  then  abruptly  switches  sign  and  decreases. 
The  large  values  of  the  derivative  in  the  boundary  condition  lead  to  a 
singular  matrix  and  the  solution  procedure  fails. 

Several  attempts  to  remedy  the  saddle  point  problem  were  made. 
Simply  neglecting  the  dependence  of  H  on  g  made  the  convergence 
procedure  much  more  difficult,  and  the  solution  would  often  diverge  only 
a  short  distance  downstream.  As  an  alternative,  an  ad  hoc  update 
procedure  for  H  was  implemented,  using  the  following  formula: 
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where  w  was  chosen  as  a  constant  less  than  unity,  generally  between  0.5 
and  0.7.  This  modified  update  procedure  was  used  when 


This  ad  hoc  update  procedure  did  limit  the  increasing  iteration 
count  as  the  solution  approached  separation.  The  procedure  does  not, 

JTJ 

however,  recognize  the  change  of  sign  in  the  derivative  —  .  The 

dp 

location  of  the  saddle  point  is  difficult  to  determine  during 
computation,  since  the  solution  cannot  identify  the  sign  switch  from  an 
upstream  position.  The  location  of  the  saddle  point  also  seems  to  vary 
with  different  flow  situations,  sometimes  immediately  at  separation  and 
sometimes  ahead  of  separation.  Since  the  location  of  the  saddle  point 
cannot  be  accurately  determined  by  the  ad  hoc  procedure,  the  method  also 
resulted  in  divergence  of  the  solution  or  the  abrupt  jump  to  another 
(incorrect)  solution  curve. 

An  additional  noteworthy  feature  of  the  hybrid  method  is  that  the 

error  between  the  specified  displacement  thickness,  6<.p,  and  the 

...  .  . 

displacement  thickness  integrated  from  the  solution,  5_„,  increases 

IN 

linearly  as  the  solution  proceeds  downstream.  This  resulted  in  a  slow 
increase  in  the  solution  error  and  the  number  of  iterations  required  for 
convergence.  It  also  lead  to  slight  differences  in  the  solutions  for 
wall  shear  methods  and  the  hybrid  method,  although  the  discrepancies 
were  acceptable. 

The  combination  of  these  problems  indicates  that  a  substantial 
investigation  of  the  hybrid  method  is  needed.  Before  the  method  can  be 
applied  constructively,  the  saddle  point  problem  must  be  solved  and  the 
method  more  fully  understood,  although  its  usefulness  is  questionable 
with  the  success  of  the  displacement  thickness  method. 
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4.1.3  Displacement  Thickness  Method 

The  displacement  thickness  method  proved  to  be  the  most  successful 
of  the  inverse  methods  applied  in  this  study.  In  addition  to  performing 
well  on  the  cases  treated  by  the  wall  shear  methods,  the  displacement 
thickness  method  obtained  very  good  results  for  cases  which  proved  to  be 
too  severe  for  the  wall  shear  methods.  In  particular,  it  was  used 
successfully  on  both  analytic  displacement  thickness  distributions 
devised  by  Carter  [13],  for  which  the  wall  shear  methods  were 
inaccurate. 

The  ability  of  the  displacement  thickness  method  to  handle  larger 
separated  regions  more  accurately  than  the  wall  shear  methods  seems  to 
stem  from  the  use  of  wall  boundary  conditions  as  discussed  previously  in 
subsection  4.1.1.  The  displacement  thickness  method  implements  a  two- 
point  spline  boundary  condition  at  the  wall  which  links  the  variables  f, 
u,  and  G  at  j  =  1  and  j  =  2.  This  additional  control  on  the  splines 
prevents  the  "wiggle"  behavior  at  the  wall  which  is  observed  in  some 
wall  shear  method  calculations. 

Several  interesting  observations  were  made  during  the  use  of  the 
displacement  thickness  method.  On  three  occasions,  a  mild  streamwise 
instability  was  noted- when  three-point  backward  differences  were  used. 
This  instability  was  first  noticeable  through  oscillations  present  in 
the  calculated  normal  component  of  velocity,  which  is  dependent  upon  f-. 
For  these  cases,  it  was  necessary  to  switch  to  two-point  backward 
differences  over  the  region  where  separation  was  present.  The  use  of 
two-point  differences  yields  accurate  solutions  through  separation 
without  difficulty. 


Another  unexpected  result  is  the  overshoot  behavior  of  the  solution 
following  reattachment  for  the  Carter  analytic  displacement  thickness 
cases.  As  discussed  in  the  previous  chapter,  the  calculations  should 
return  to  a  flow  situation  with  a  near  zero  pressure  gradient,  as  for  a 
flat  plate.  For  the  present  method,  however,  the  calculations  produce  a 
flow  situation  which  continues  to  accelerate  following  reattachment, 
overshooting  the  flat  plate  value  of  6  i/Re.  The  overshoot  problem, 
while  present  for  both  displacement  thickness  distributions,  is  much 
more  pronounced  for  the  more  severe  distribution. 

The  suspected  cause  for  this  behavior  is  streamwise  discretization 
error,  which  is  aggravated  by  the  use  of  two-point  backward  differences 
for  the  more  severe  displacement  thickness  case.  As  a  check,  the 
streamwise  steps ize  was  halved  and  the  cases  rerun,  as  discussed  in  the 
previous  chapter.  The  calculations  with  halved  streamwise  stepsize  show 
a  marked  improvement  over  the  initial  calculations,  although  the 
overshoot  behavior  is  still  noticeable.  The  method  has  proven  to  be 
sensitive  to  strongly  accelerating  flows,  however,  and  this  sensitivity 
may  be  the  cause  of  the  remaining  discrepancies. 

The  present  formulation  was  also  found  to  be  very  sensitive  to 
input  data  such  as  displacement  thickness  or  wall  shear  distributions. 
This  is  clear  after  examining  case  number  6,  treating  the  Horton 
displacement  thickness  distribution  obtained  experimentally  by  Ntin 
[11].  Here,  the  graphical  edge  velocity  distribution  given  by  Horton 
was  fitted  using  cubic  spline  interpolation,  and  a  direct  solution  was 
computed  to  obtain  the  displacement  thickness  distribution  to  a  location 
near  the  separation  point.  We  note  that  the  edge  velocity  distribution 
is  matched  at  two  locations,  x  ■  0.58  and  x  =  0.62.  This  is  necessary 


due  to  the  graphical  presentation  of  the  edge  velocity  in  Reference 

[111. 

Once  the  displacement  thickness  distribution  is  obtained  from  the 
direct  solution,  it  is  fitted,  again  using  cubic  spline  interpolation, 
to  the  distribution  given  by  Horton  for  the  separated  region  of  the 
boundary- layer  flow.  The  matching  procedure  is  completed  over  an 
interval  from  x  =  0.61  to  x  =  0.66.  A  portion  of  the  direct 
displacement  thickness  distribution  is  also  smoothed  between  x  =  0.57 
and  x  =  0.585.  Although  care  is  taken  to  produce  continuous 
distributions  of  ug  and  <5v/Re,  the  derivatives  of  these  quantitied  are 
not  forced  to  be  continuous  at  the  matching  location. 

The  inverse  solution  procedure  proved  to  be  very  sensitive  to  the 
derivative  discontinuities,  as  is  shown  in  Figure  38.  The  locations  of 
the  matching  points  can  be  clearly  located  by  the  peaks  in  the  skin 
friction  coefficient  distribution.  From  this  result,  it  is  clear  that 
for  the  present  formulation,  great  care  must  be  taken  to  assure  that  the 
input  data  is  as  smooth  and  free  of  numerical  noise  as  possible  to  allow 
for  an  accurate  calculation.  The  situation  here  is  aggravated  by  the 
lack  of  available  tabular  data. 

4.2  Direct  Method 

The  direct  solution  option  of  the  present  boundary- layer 
formulation  was  tested  in  cases  1,  2,  and  3,  and  used  to  obtain  a 
displacement  thickness  distribution  in  case  6.  It  performed  well  for 
all  laminar  test  cases  considered  here,  with  quadratic  convergence 

f  I 

consistently  observed.  Upon  reaching  negative  values  of  f^  ,  the 
solution  diverged  as  expected  for  standard  direct  boundary- layer 


solutions . 


The  present  direct  method  is  interesting  since  only  two  boundary 
conditions  must  be  changed  to  switch  the  solution  procedure  from  inverse 
to  direct  mode.  To  obtain  a  direct  solution  with  the  present 
formulation,  the  additional  condition  on  (3  must  be  enforced  and  the 
inverse  condition  discarded.  The  direct  method  retains  the  4X4  block 
structure  of  the  inverse  methods.  The  mechul  function  equation  on  3 
becomes  a  "dummy"  equation  since  3  is  specified  and  the  equation  is 
satisfied  before  iteration  begins.  This  suggests  that  the  solution 
procedure  could  be  easily  switched  from  direct  to  inverse  mode  during 
calculation,  although  this  was  not  attempted  during  this  study. 

With  the  edge  velocity  distribution  specified,  3  is  obtained  using 
finite  difference  formulas  obtained  from  the  Lagrange  interpolation 
polynomial.  In  general,  three-  or  four-point  backward  differences  are 
used.  Near  the  starting  solution,  forward  differences  are  used  if  the 
flow  is  laminar  several  stations  downstream.  If  the  boundary  layer 
becomes  turbulent,  3  must  be  updated  after  each  iteration,  since  it  is 
dependent  upon  the  edge  value  of  pt. 

Additional  comments  on  the  direct  method  are  included  in  the 
following  sections  4.3,  Boundary  Conditions,  and  4.4,  Turbulence 
Application. 

4 . 3  Boundary  Conditions 

Throughout  the  present  investigation,  the  methods  developed  have 
proven  to  be  very  sensitive  to  the  boundary  conditions.  A  prime  example 
is  the  results  obtained  by  the  wall  shear  method  for  certain  separated 
regions,  as  discussed  in  section  4.1.  This  behavior  was  completely 
unexpected,  since  the  boundary  conditions  applied  at  the  wall  are  quite 
typical.  The  use  of  the  splines,  however,  is  not  typical. 


Unexpected  results  related  to  the  two-point  spline  relations  were 
observed  in  other  situations  aside  from  the  behavior  of  the  wall  shear 
methods  in  separated  flows.  The  oscillations  produced  by  transferring 
the  two-point  condition  on  3  from  the  wall  to  the  far  field  are  another 
instance  of  such  unexpected  behavior.  Although  the  boundary  condition 
was  expected  to  be  satisfactory  at  the  far  field,  the  resulting 
oscillations  caused  the  solution  to  diverge  at  a  streanrwise  location 
short  of  the  distance  obtained  when  the  condition  on  3  is  implemented  at 
the  wall. 

Similar  behavior  was  observed  for  the  direct  method.  Replacing 
Eq.  (66a)  with  Eq.  (66c),  both  of  which  should  perform  satisfactorily  in 
the  far  field,  caused  the  u-velocity  profile  near  Hoc  to  overshoot  the 
specified  value  of  unity  at  Ho,,.  This  overshoot  tended  to  increase  as 
the  solution  proceeded  downstream,  and  caused  a  decrease  in  solution 
accuracy. 

Other  instances  of  such  behavior  were  observed  in  the  current 
investigation,  but  only  the  most  severe  cases  are  discussed  here.  In 
all  instances,  however,  the  solution  procedure  was  found  to  be  extremely 
sensitive  to  the  particular  implementation  of  the  two-point  spline 
relations.  Often  it  was  found  that  one  condition  adversely  affected  the 
solution  while  another  provided  accurate  results,  although  both  were 
expected  to  be  adequate.  The  cause  of  this  sensitivity  is  unclear  from 
the  results  of  the  current  investigation. 

4.4  Turbulent  Application 

As  noted  in  the  previous  chapter,  the  attempt  at  obtaining  a 
solution  for  a  turbulent  boundary  layer  was  far  less  successful  than  the 
laminar  solutions.  As  a  test  case,  the  present  direct  and  inverse 
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methods  were  used  to  compute  a  turbulent  boundary  layer  on  a  flat  plate. 
As  the  solution  proceeds  downstream,  however,  the  accuracy  in  the  far 
field  rapidly  deteriorates.  The  value  of  the  shear  ( f ' ' )  near  the  edge 
grows  steadily  with  downstream  distance,  while  the  values  of  u  ( f ' )  near 
Hoo  overshoot  the  edge  value  of  unity.  Oscillations  in  the  far  field 
follow,  leading  to  divergence  of  the  solution  only  a  short  distance 
downstream. 

The  far  field  boundary  conditions  were  initially  suspected  as  the 
cause  of  this  phenomena,  since  the  methods  had  previously  been  shown  to 
be  sensitive  to  the  particular  two-point  spline  conditions  implemented. 
While  changing  the  particular  implementation  of  the  two-point  spline 
conditions  produced  minor  adjustments  in  the  solution,  such  as  slightly 
smaller  values  of  shear  near  the  edge,  the  modifications  did  not  prevent 
the  oscillations  or  the  early  divergence  of  the  solution. 

At  present,  two  possibilities  exist  which  may  be  the  cause  of  the 
turbulent  solution's  behavior.  The  most  likely  cause  is  the  presence  of 
a  streamwise  instability  in  the  present  formulation  which  has  a  minimal 
effect  on  laminar  calculations,  but  which  is  much  more  important  for 
turbulent  calculations.  This  is  a  prime  suspect  since  streamwise 
fluctuations  were  found  in  the  turbulent  solution  attempts. 

Oscillations  were  also  noted  in  the  normal  velocity  profile,  similar  to 
those  found  for  the  laminar  displacement  thickness  cases.  For  the 
turbulent  case,  switching  to  two-point  backward  differences  prolongs 
solution  convergence  by  several  streamwise  stations,  but  does  not 
prevent  the  solution  from  diverging. 

The  other  possibility  is  that  the  present  formulation  is  not 
compatible  with  the  use  of  the  modified  Levy-Lees  transformation. 


Although  this  possibility  is  more  remote,  it  cannot  be  rejected  at  this 
stage.  It  is  also  possible  that  the  behavior  of  the  turbulent  solution 
is  the  result  of  a  combination  of  these  causes.  The  actual  cause  cannot 


be  definitely  identified  without  further  investigation. 

Because  a  turbulent  solution  was  never  successfully  completed,  the 
adaptive  grid  coordinate  in  its  present  form  was  never  fully  tested. 

The  laminar  section  of  the  coordinate  was  tested  on  several  laminar 
cases,  yielding  results  typical  of  a  hyperbolic  normal  distribution. 
Numerical  results  for  these  laminar  cases  are  not  included  here  since 
they  are  not  noteworthy. 


A . 5  Advantages  and  Disadvantages  of  Similarity  Transformations 

Several  advantages  of  the  modified  Levy-Lees  transformation  have 
been  discussed  in  Chapter  Two,  such  as  the  ability  to  capture  the 
laminar  and  turbulent  boundary- layer  growth  and  to  provide  a  self¬ 
starting  solution.  This  transformation,  or  the  more  familiar  original 
Levy-Lees  form,  has  been  a  standard  part  of  many  direct  boundary- layer 
solutions.  For  inverse  flows,  the  use  of  the  transformation  is  not 
standard.  In  most  cases,  the  streamwise  coordinate  x  is  left  alone, 
although  the  normal  coordinate  may  be  put  into  similarity  form.  This  is 
mainly  because  the  transformed  streamwise  coordinate  E,  is  a  function  of 
the  edge  velocity,  which  is  obtained  during  an  inverse  solution,  rather 
than  known  beforehand  for  a  direct  solution. 

Leaving  the  streamwise  coordinate  untransformed,  however,  precludes 
the  use  of  the  modified  Levy-Lees  transformations  and  its  desirable 
effect  of  capturing  the  turbulent  boundary- layer  growth.  This  study  has 
shown  that  with  proper  linearization,  use  of  the  transformed  streamwise 
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coordinate  need  not  result  in  additional  iterations,  particularly  for 


laminar  flows. 

For  direct  laminar  solutions,  the  use  of  the  modified  Levy-Lees 

transformation  is  straightforward,  since  ue  is  specified  and  pfc  reduces 

e 

to  the  molecular  viscosity.  For  direct  turbulent  solutions,  the  use  of 

the  transformation  becomes  more  complicated  since  p^_  must  be  estimated 

e 

and  updated  at  each  streamwise  station  because  of  its  dependence  on  the 
eddy  viscosity,  e. 

Inverse  solutions  are  complicated  further  since  ue  is  obtained 

during  the  solution  procedure.  Further,  for  inverse  turbulent 

solutions,  5  must  be  updated  with  u  and  p^.  as  both  change  during  the 

u  e 

solution  iteration.  Although  this  additional  update  procedure  may  tend 
to  slow  the  convergence  rate,  this  inconvenience  is  outweighed  by  the 
advantages  provided  by  applying  the  transformation.  Additional 
iterations  required  for  convergence  can  also  be  reduced  by  proper 
linearization,  as  discussed  in  section  2.9.  Unfortunately,  due  to  the 
lack  of  success  with  the  turbulent  solution  in  this  study,  the 
usefulness  of  the  modified  Levy-Lees  transformation  could  not  be 
verified  beyond  its  simplification  of  equation  form  and  algorithm 
development. 

4 . 6  Pivoting 

The  present  direct  and  inverse  methods  employ  partial  pivoting  in 
the  L-U  decomposition  to  prevent  the  buildup  of  roundoff  error.  In  an 
effort  to  determine  whether  an  ideal  ordering  of  the  equations  exists, 
the  pivoting  indices  in  the  block  matrix  solver  were  observed  during 
|  solution  of  a  nonsimilar  boundary- layer  flow.  By  checking  these 

indices,  it  is  possible  to  determine  the  pivoting  elements  of  the  L-U 

i 
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decomposition  procedure.  If  the  decomposition  routine  pivots 
consistently  about  the  same  elements,  an  ideal  order  of  the  equations  in 
block  matrix  form  is  implied. 

The  check  of  the  pivoting  indices,  howe\  r,  shows  no  such 
consistencies.  Not  only  does  the  pivoting  order  change  continually 
across  the  boundary  layer,  but  it  also  changes  as  the  solution  proceeds 
downstream.  This  indicates  tnat  no  predetermined  order  of  the  equations 
exists.  Thus,  regardless  of  the  equation  ordering,  pivoting  is  an 
essential  process  for  obtaining  an  accurate  solution. 


Chapter  5 


CONCLUSIONS  AND  RECOMMENDATIONS 

A  unified  approach  to  calculating  direct  and  inverse  solutions  for 
two-dimensional  boundary  layers  has  been  presented.  By  a  simple 
interchange  of  boundary  conditions,  the  method  is  capable  of  providing  a 
direct  solution,  with  a  pressure  distribution  specified,  or  one  of 
several  inverse  solution  options,  depending  upon  the  specification  of 
wall  shear,  skin  friction  coefficient,  or  displacement  thickness 
distributions.  A  modified  form  of  the  mechul  function  method  is  used  to 
obtain  the  inverse  solutions.  Spline/finite  difference  discretization 
yields  a  fully  implicit  scheme  which  converges  quadratically  for  laminar 
flows,  when  proper  linearization  is  applied. 

For  laminar  flows,  the  method  works  quite  well,  although  it  was 
discovered  that  the  wall  shear/ skin  friction  coefficient  options  do  not 
converge  properly  for  separation  regions  with  backflow  velocities 
greater  than  one  percent  of  the  edge  velocity.  This  behavior  is  due  to 
the  nature  of  the  splines  combined  with  the  lack  of  a  sufficiently 
strong  two-point  boundary  condition  at  the  wall.  Without  proper  control 
through  the  wall  conditions,  the  splines  allow  oscillations  in  the 
solution  profiles  in  regions  of  separation.  The  displacement  thickness 
option,  which  does  employ  such  a  wall  boundary  condition,  is  capable  of 
treating  regions  of  separation  with  larger  backflow  velocities 
(approximately  ten  percent  of  the  edge  velocity). 

The  solution  procedure  is  also  sensitive  to  the  use  of  the  three- 
point  backward  differences,  sometimes  requiring  the  use  of  two-point 
backward  differences  for  accurate  results.  For  laminar  flows,  this 
sensitivity  appeared  during  the  calculation  of  separated  regions,  and 
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was  manifested  by  oscillations  in  the  normal  component  of  velocity.  It 
is  believed  that  this  sensitivity  is  connected  to  the  instability  noted 
during  the  attempted  turbulent  calculations.  In  addition  to  the  normal 
velocity  oscillations,  streamwise  fluctuations  were  also  noted  in  the 
turbulent  cases.  This  instability  resulted  in  the  failure  of  the 
turbulent  calculations.  It  is  unclear  why  this  instability  strongly 
affects  the  turbulent  case  while  having  little  effect  on  the  laminar 
results . 

These  conclusions  cast  some  doubt  on  the  utility  of  spline 
discretization  for  all  but  a  small  portion  of  boundary- layer  flows,  and 
particularly  for  cases  concerning  regions  of  separation  or  turbulence. 
The  sensitivity  of  the  present  method  to  the  two-point  spline  boundary 
conditions  contributes  strongly  to  this  conclusion.  It  should  be  noted, 
however,  that  the  spline  formulation  performed  excellently  in  many 
instances,  yielding  very  accurate  results  with  a  minimum  of  mesh  points. 
The  use  of  a  strong,  two-point  boundary  condition  at  the  wall  is 
essential  for  solutions  involving  reversed  flow.  Satisfactory  results 
have  been  obtained  only  when  such  a  condition  is  used.  Unfortunately, 
this  form  of  boundary  condition  is  sometimes  difficult  to  implement,  as 
demonstrated  by  the  wall  shear/skin  friction  options. 

The  work  of  Rubin  and  Khosla  [44]  indicates  that  turbulent,  non¬ 
separating  boundary  layers  can  be  successfully  treated  by  spline 
discretization.  Therefore,  prior  to  completely  discarding  the 
spline/finite  difference  formulation,  a  full  stability  analysis  must  be 
performed  to  determine  the  source  of  the  streamwise  instability  which 
prevents  the  calculation  of  accurate  turbulent  results  and  impedes  the 
calculation  of  some  separated  solutions  for  this  method  formulation. 


2 


Farther  research  into  the  application  of  an  adequate  boundary  condition 
for  the  wall  shear  option  is  also  required  to  establish  the  minimum 
requirements  of  such  a  condition  and  how  it  could  best  be  implemented. 


Sr  ^  V*. 
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Appendix  A:  BACKWARD  DIFFERENCES  WITH  VARIABLE  STEPSIZE 
Streamwise  differencing  with  variable  stepsize  is  essential  for  the 
present  inverse  method  since  the  transformed  streamwise  coordinate,  £  is 
determined  with  the  solution.  Thus,  although  a  constant  stepsize  in  x 
may  be  specified  with  the  input  data,  the  values  of  A?  between  each 
streamwise  station  will  vary.  To  obtain  accurate  difference  formulas 
with  variable  stepsize,  the  Lagrange  interpolation  formula  is  used  to 
derive  three-  and  two-point  backward  differences. 

The  quadratic  Lagrange  interpolation  formula  is 


(x-x  )(x-x_)  (x-x  )(x-x_) 

f(x)  =  f(x.)  7 - 77 - 7  +  f(x_)  7 - 77 - r 

1  (x^x^Mx^-x^)  2  (x2*x1)(x2-x3) 


(x-x. )(x~x9) 

+  f(x~)  t - w — - — 7 

3  (x^-Xj )(x--x2) 


(A.  1) 


Letting  f(x)  *  y  -  y^,  Eq.  ( A . 1 )  becomes 


(A. 2) 


(A. 3) 


where 
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A1  =  (x-x^)  +  (x-x^) 


=  (x-x^)  +  (x-x2) 


In  terns  of  the  streamwise  index  i,  the  three  point  backward  difference 
formula  at  x  is 


(S'i  ■  57  (yi-i  ■  +  5;  <yi-2 '  V 


(A.  A) 


where 


A.  =  x.-x.  - 

1  ii-2 


A_  =  x. -x.  , 
2  11-I 


Bi  ■  (*i-rxi)(vrxi-2) 

b2  -  (xi.2-xi,(xi.2-xi.1) 

Equation  (A. 4)  can  be  written  in  the  alternative  form 


,dv,  A1  A2.  A1  A2 

"  ‘  (B[  +  B?  h  +  B[  yi-l  +  yi-: 


(A. 5) 


(-p) .  =  ay.  +  by.  .  +  cy. 
dx  i  J 1  •'i-l  J 1- 


(A.6) 


Equation  A. 6  is  the  form  used  to  discretize  the  streamwise 


derivatives . 


For  a  two-point  backward  difference  formula,  Eq.  (A. 4)  reduces  to 


rdv-,  ”1  ,  » 

(S'i  ‘  (yi-l  '  yi> 


(A. 7) 


where 


A1  "  1 


B,  ■  x.  , -x. 
1  i-l  1 
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Appendix  B:  TRANSFORMED  EDDY  VISCOSITY  FORMULAS 
The  untransformed  eddy-viscosity  model  of  Baldwin  and  Lomax  [24] 
described  in  Section  2.6.  For  this  problem,  the  eddy-viscosity 
equations  must  be  non-dimensionalized  and  transformed  to  (L,  n) 
coordinates.  The  following  is  a  complete  list  of  all  transformed 
quantities  related  to  the  turbulence  model. 

Inner  region: 


where 


5.  =  /Re  k2  /2£  q2  l1  jf  |  v 
i  c  1  nn  tr 


l  =  1  -  exp  (- 
C  A 


y+  =  RelM  n  [/2l  jf  j  ] 

J  1  nn  w 


1/2 


Outer  region: 


Eo  Re  K  CCP  FWAKE  FKLEB  (n)  Ytr 


where 


F  =  < 
WAKE  \ 


Re"1^2  n  F 

u  'MAX  max 


D  -1/2  „  /2?  2  , 

Re  CWK  u  nMAX  UDIF/FMAX 
e 


+ 


F(n)  =  u  n J f  J  [1  -  exp  (-  ■*-)] 
e  nn  A  + 


Fu,v  =  maximum  [F(n)l 
MAX 


(B.l) 

(B.  2) 

(B.  3 ) 

(B.  4) 

\ 

(B.5) 
\  the  smaller 

(B.6) 

/ 

(B.  7) 


'MAX 


MAX 


is 


-  n  where  F(n)  =  F, 
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Appendix  C:  SPLINE  AND  BLOCK  TRIDIAGONAL  COEFFICIENTS 


In  block  matrix  form,  the  system  of  spline  equations  for  the  nodal 


points  2  <  j  <  N  at  a  streamwise  station  i  are  as  follows: 


where 


B.  Z.  .  +  A.  Z.  +  C.  Z.,,  =  R. 
J  J-l  J  J  J  J+l  J 


-r  a. 


-  B 


°2  * 

—  bj-i(ej-rI) 
j-i 


1  *, 

—  b  (e  -1) 

N.  3  J 
J 


Nj Gj 


TO 
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cr  =  u —  >  h.  =  N.  -  N.  ,  or  h .  =  n .  -  n  .  , 

hj  J  J  J-l  J  J  'j-l 


A .=  C_.  b.G. 
J  2j  j  j 


E.  =  2u .  (p.  +  0.  C,  )  +  R  .  ( 0 . - 1 )  b?  G. 
J  J  J  J  J  J  J  J 


C.  -  bj(i-ej)[C2  f.  +  Rj(0j-l)Uj  +  C3  ] 


D.  =  C.  -  r1  — 
J  J  hj  °2 


e.  -  c.  -  r1  — 

j  i  ijij 


F.  =  c-xi  -  IT 

J  J+l  hj  5 


G.  =  uT  -  1 
J  J 

ne 

R.  =  2  5(^). 

J  Nn  j 

The  backward  differencing  coefficients  are  defined  as 
Cl..  =a5i 
C2  =  *(2a£.  +  1) 


C3 .  .  25i  <bfi-l  + 

1J  J 


V,  ■  ?i  <bui-l  +  cui'2)j 

1J  J 

with  a,  b,  and  c  as  given  in  Appendix  A.  The  eddy  viscosity 


.v.»r 


»’,*L 

i*. 


coefficients  are 


*  -1 
bj  =bj 


b.  =  — 

J 


ut  =  1  +  e 


u  =  1  +  e 
Mt  o 

e 


2z .  +  1 
i 


for  inner  region 


for  outer  region 
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